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HEAT TRANSFER BY LAMINAR FREE CONVECTION 
IN ENCLOSED PLANE GAS LAYERS 


By G. POOTS (Department of Theoretical Mechanics, University of Bristol) 


[Received 31 May 1957] 


SUMMARY 

A numerical method based on the use of orthogonal polynomials is given for the 
solution of two simultaneous partial differential equations which govern the heat 
transfer by laminar free convection in enclosed horizontal, oblique, and vertical gas 
layers. The heat transfer through the enclosed gas layer depends only on the Rayleigh 
nuinber A = g(T,— T,)d*/(T, kv), the Prandtl number o = v/k, the aspect ratio of the 
rectangular cell a = I/d, the angle of inclination ¢ of the hot and cold walls to the 
vertical, and finally the temperature distribution in the remaining two walls or 
border strips. 

Preliminary thermal results are given for air taking o = 0-73, a = 1, ¢6 = 0 and 
A = 5x 10%, 10%, 2-5 10%, 5 10°, and 10‘, and the temperature distribution along 
the border strips to be T = T)+(7,—%)y/d. A comparison of the theoretical pre- 
diction of the heat transfer across the cell with the experimental data of Mull and 
Reiher (1), as correlated by Jakob (2), can only be made in the region A = 104, 
Calculated results and empirical formulae obtained by Jakob agree favourably in 
this region considering that all the measurements of Mull and Reiher were made for 
large values of the aspect ratio and virtually non-conducting border strips. 

Streamlines and isothermals have been plotted for A = 10* which clearly indicate 
the existence of an isothermal core in the cavity having constant vorticity, as postu- 
lated by Batchelor (3). 


1. Introduction and basic theory 


Tue heat transmission through horizontal, oblique, and vertical enclosed 
gas layers was investigated experimentally by Mull and Reiher (1) and the 
thermal results correlated by Jakob (2). Batchelor (3) derived the equations 
governing heat transmission across an air space between two vertical plane 
boundaries, distant d apart, and held at different temperatures, the air 
space being closed by two horizontal plane boundaries, or border strips, 
distant / apart. Quantitative information was given on the rate of heat 
transfer across the cold boundary for large aspect ratios « = I/d. It is the 
purpose of the present paper to give a numerical method for solving the 
two simultaneous partial differential equations which govern the heat 
transmission in a closed rectangular cell having hot and cold walls inclined 
at an angle ¢ to the vertical. The boundary condition on the border strips 
is assumed to be either linear and of the form 7’ = 7,+(7,—7)y/d, or that 
there is no heat loss into the surrounding medium, namely @7'/éx = 0, 
Consider an air space enclosed between two plane parallel boundaries 


{ Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 3, 1958) 
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inclined at an angle ¢ to the vertical, and by two border strips of length d 
and distance / apart as shown in Fig. 1. Let 7, and 7, be the absolute 

temperatures of the hot and cold walls respectively. 
y Assuming that 





(i) temperature difference 7,—7, is small com- 
pared with the absolute temperature 7, of 
the cold boundary, 

(ii) viscosity, density, and thermal conductivity 
are constant except for the effect of the den- 
sity variation in producing buoyancy force, 

d \ (iii) fluid is incompressible and viscous heat dissi- 
pation may be neglected, 

(iv) motion is two-dimensional, 


Ty 


_—Te 














Fic. 1 


then subject to the above restrictions the equations governing free con- 
vection by laminar flow in the cell may be formulated (Goldstein (4)). Thus 
if the coordinates x, y are chosen as in Fig. 1, and the corresponding velocity 
components are w, v, then the basic equations of motion are 





ou eu 1 Op T—T, 

Uu— ) = —— Rs. KS 2 
ont Oy =~ ie geos T, -vV2u, (1.1) 
ov ov 1 Op . ,(T—T, 

u— ae GER enw ae ae ee Soeinbtecabel V2) 9 
ee * ay Po &Y asin g| qT) aaa . wing 
cu ov oe 

ta By == 0, (1.3) 


ns. 
if lated eae Zs (1.4) 


where the suffix 0 signifies conditions at the cold wall, p is the fluid densit y; 
p the pressure, g the gravitational constant, v the kinematic viscosity, and 
k the thermometric conductivity. The boundary conditions to be imposed 
on the variables wu, v, and 7' are 


u=v=0, T=T% aty=0 for 0 <2 <l, 
we=eo=0, T=7T, aty=d for 0 <2 <l, 
u=v=0 atz=0,l for0 <y <d, 

T = T%+(T,—T)y/d or @T/éx=0 atx = 0,1 for0 <y <d. 


Following Batchelor (3) by taking d as a representative length for the 
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coordinates x and y (without changing notation), and defining a dimension- 
less temperature distribution and stream function by 


6=(T—T,)(T.—T,), w= ka, s—-eee aN 


ea 
respectively, then on eliminating the pressure p the governing equations 
become 


V9 = — T_— =, (1.6) 


v4, = A(cosg 5 —sing =) + 2/0) ! (vey) (1.7) 


o\ex cy ,’ 
where o is the Prandtl number, and A = oG = (T,—7,)gd*/(T, kv) are the 
Rayleigh and Grashof numbers respectively. The associated boundary 
conditions are 


mye ies = 
at he 6=0 aty=90, 


yp 


cy 


b= 0, 6=1 aty=1 for0 <2 <a, 


y= ~ == 0, 6 = y (Case I) or <- 0 (Case IT) atxr=0,a 
Ce 


forO <y <1. 
The form of the governing equations (1.6) and (1.7) and the above boundary 
conditions shows that the quantities a, ¢, A, and o, together with the border 
strip boundary conditions, are sufficient to determine % and @ uniquely. 
The thermal quantity measured experimentally is the rate of heat trans- 
mission through the cold boundary. If this rate is Q heat units per second 
per unit depth of boundary in the z-direction, then a dimensionless quantity 
describing the heat transfer is the Nusselt number 


N = ae = a d . 1.8 
k(T,—T) . CY} y=0 z ie 


As has been shown by Batchelor (3) approximations to both 6 and % may 
be made enabling quantitative information on N to be obtained for the 
various types of flow occurring in vertical cavities (6 = 0). For A < 108 
Batchelor develops a power series expansion in terms of A for @ and % as 
follows 


b= > AM(Cy), O=y+ > A°(2,y). (1.9) 


On substitution of (1.9) into the governing equations (1.6) and (1.7) and 
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equating coefficients of like powers of A, he reduces the problem to the 
solution of a series of linear partial differential equations in #,,(x,y) and 
6,,(x,y). Due to the symmetry properties of the @,,, 


x 
\ 


a * (00, 
N= a+ > Yan _ ms — | ) az (1.10) 
—_ | J\ ay }y-0 | 
0 

The value of y, was estimated to be of the order 10-* when a = I/d was 
not too different from unity. A section is then given on the case of a 
general value of A as «> 00. Here Batchelor argues that in the region not 
near the ends of the cavity @ and % take up their asymptotic forms 


6 = y, y= Sux) (1.11) 


and since in this region the rising stream loses heat at a rate equal to the 
net flux of heat across the horizontal plane, he then obtains 

N~ + ("55 4 (1.12) 

d‘ \ 720 

where f is an undetermined constant taken to be unity. In a final section 
on the general value of « as A - 0c Batchelor postulates that inside the 
cavity there exists an isothermal core having temperature @ = 4, and having 
constant vorticity. Boundary layer equations are set up for the surrounding 
continuous boundary layer in which temperature and velocity make rapid 
transitions to their assigned values at the boundary. Batchelor found that 
since # and the velocity components oscillate in the boundary layer these 
quantities could not be satisfactorily represented by polynomials of small 
degree, nor did the equations yield to linearization of the Oseen type. 
However, although not solving these equations, Batchelor has obtained 
estimates of the rate of heat transfer by assuming that the flow near each 
wall bears resemblance to either forced convection of a stream or free 
convection past a plate under certain boundary conditions. 

In section 2 of this paper a technique, based on the use of orthogonal 
polynomials, is developed for obtaining numerical approximations to the 
solution of the two partial differential equations (1.6) and (1.7). Numerical 
solutions are also tabulated, which have been obtained on application of 
the technique to the case of a cavity having « = 1,46 = 0,andA = 5X 102, 
108, 2-5 108, 5 108, and 104. Details of the calculations involved are 
given in Appendixes A and B. From these solutions thermal results may 
be extracted which are given and compared with existing experimental data 
in section 3. Finally, in section 4 a statement of conclusions is made on the 
numerical technique used, and on the thermal results so far obtained. 
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2. Development and application of the numerical technique 

The purpose of this section is to show how, by the use of orthogonal 
polynomials, the solution of the governing equations (1.6) and (1.7) may 
be reduced to the numerical solution of two sets of coupled algebraic 
equations. For convenience in notation we convert the region of integration 
from that of a rectangle x = 0, « and y = 0, 1 to a square region of unit 
dimension by the substitutions 


wt=af y=; (2.1) 
also to standardize the boundary conditions on @ we introduce 
HE, n) = n+O(€, 0). (2.2) 


The governing equations (1.6) and (1.7) then become 


, 120 20 _ 1/0 2% ap 
V20 = -— —_ =e i, 
a? og a ! én On C& G€ 





(2.3) 


. a ao) 1 {2 0. a S. v2 Oyp\_ 
Viy = sous’. —+cos¢ “sing? )+ + colae V8 xb) = 4 a¥) 2" 
The implied boundary conditions on ¢ are (2.4) 
$=7 = , atyn=—90,1 forO <é<l, 

“ 9 
oars até=0,1 forO0<y <1. (2.5) 
Cc 


Furthermore, if we take the temperature distribution on the border strips 
to be linear, then 


O=0 atyn=0,1 fordO<é< 1, 
and at €=0,1 for0 <7 <1 (Case I); (2.6) 
or if the border strips are perfect insulators, then 

O=0 atyn=0,1 ford <é <1, 


and —o =0 at&é=0,1 for0 <7 <1 (Casell). (2.7) 


Case I. Linear temperature distribution on the border strips 


We assume that © and Y may be represented by the complete double 
series of orthogonal functions 


8= 4 4 A, _Sin pr€é sin q7 (2.8) 
and b= > Brg Xpl€Xaln), (2.9) 
o=1 


where the A,, and B,,, are constants to be evaluated. 
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The double Fourier series expansion (2.8) satisfies the boundary condi- 
tions (2.6); also the expansion (2.9) satisfies the required boundary con- 
ditions (2.5), provided the functions X,,(£) are chosen to satisfy the fourth 
order Sturm—Liouville equation 


d‘X,, . ; 
dé ta up X, (2.10) 
with X (0) = X,(1) = X},(0) = X,(1) = 9, 


and where a prime denotes differentiation with respect to €. The differen- 
tial equation (2.10) is the well-known equation determining the modes of 
vibration of a uniform beam clamped at both ends. The appropriate eigen- 
functions may be given in the form 
X,(€) = (cosh p, € — cos p, €)—8,(sinh nw, € — sin p, €), (2.11) 
if the eigenvalues p,, satisfy the transcendental equation 
cosh 2, Co8 z,—1 = 0, 

and B, = (cosh up, — cos p,)/(sinh p, — sin p,).T (2.12) 
Furthermore the functions X,, are orthogonal in the range of integration 
(0,1), thus : 

| X,(€)X,(£) dé = 8, (2.13) 

0 


where 5,,,, is the Kronecker delta. 


We now seek the solution of (2.3) in the assumed form (2.8) where the 
Fourier coefficients are 


11 
4, = 4 | | O(é, n)sin mz€ sin nan dédy (2.14) 
0 0 


for m,n = 1, 2, 3,..., ete. On integration by parts and using the imposed 
boundary conditions on 0, we have 





1 
4 rf eo a0) 
A,.,. = | “~ cos mné sin nan dédn = id iE - sin mr€ cos nan dédy 
Fe on 


m 7 nr 
00 


g a4 
1 
~ - an) 8 Tes 


00 
0 
1 
rae Vie 


+ The values of 4», and 8, have been tabulated by Young (5), also the integral 





- oF 





(2.15) 


1 
f Xe) XUE) dé. 


The functions X,(£) for p = 1(1)6 at € = 0-0(0-05)0-5 have been tabulated by Inglis (6) 
and Corlett (7). 
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It then follows from (2.15) that on the left-hand side of the governing 
equation (2.3) 


V29 = — > 24 ni(P +0)4,, _sin pr€é sin gz, (2.16) 
p=1q ’ 
and on the right-hand side 


(o2- CO yp _ 4 
alcE én On CE OF 


_. — op -. > m A, (p= cos pr€ sin g77 = qm sin pr€ cos qn =) . 
al o&€ > i en 0& 

(2.17) 

In the usual manner if the double Fourier sine series (2.16) is equal to 
(2.17), then the Fourier transforms of the two sides must be equal. Hence 
by multiplying both sides of the equation by 4sin m7ésin ny and inte- 
grating over the region = 0, 1 and » = 0, 1 we obtain an infinite set of 
algebraic equations in A,,,,. Thus for m, n = 1, 2, 3... 


1 
(= +n ‘\A —_ == ie —sin m7€é sin nary dédn 
00 


-£5 D3 Ana {| (» cos prt singe — 


7 p=1 = 
—qsin pr€ cos qrn )sin mr€sinnan dgdn. (2.18) 


It remains now to find suitable expressions for the partial derivatives &j/e& 
and @b/én. Suppose we expand é//0€ as follows : 


B= SS ban Xml OXa(n) (2.19) 
11 
then ban = | | EX W(X qn) ded. 


On integration by parts a using the boundary conditions (2.5) we find 
11 


ban = — f \¥ PXin(E)X n(n) dédy 


= 


™" 1 


cae . 3B va | Xp(€)Xin(€) aE j X,(n)X,() dy 


By» j X,(6)X%(€) dé = p’ B,, | X(E)Xq(E) dE. 
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Thus on substituting in (2.19) we have 
1 


Fe 2 By 2, Bom | Xo(OXmE) de) Xu(€OXuln) 


of m=1n=1‘p=1 


which on pogeeeg gives 


> > B, n)( > X,,(€) j X,(€)Xn(€) dé) 
eo p=1n=1 
=> 4 By» X;(€)X (0). (2.20) 


Hence term by term aitisenntintion of (2.9) to obtain ép/a and éy/én is 
justified. Finally substituting (2. ry in (2.18) we find 


(Zmt+nt) dna = - » = i, > B,,.(a8,07,)-+Ymns (2.21) 
where wiih 
mn — mp SS 5 > Anal 2 B, APR, p,m)T(8,q,n)— 
p=1q= r=l1 s= 


—qR(s,q,n)T(r, p,m)}]. | (2.22) 
Definitions and formulae for evaluation of the integrals a’, b7,, R(r, p,m) 
and 7'(s,q,n) are given in Appendix A. 

By a similar procedure to that outlined in obtaining the result (2.20) from 
(2.9), it can be verified that the remaining partial derivatives occurring in 
equation (2.4) may be obtained by term-by-term differentiation of the 
series (2.9). On making use of the properties of the basic functions X,, the 
partial differential equation (2.4), defining %, may be reduced to the numeri- 
cal solution of a set of algebraic equations in B,, ,, by the technique employed 
in treatment of (2.3). Thus for m, n = 1, 2, 3...., 


1 Rs 
(ss Mn ¥ oe Dann Dyn +u8) an 


= AC, ,co8g—45 > Ap, ap bp cos —* az bz sing) — 


p=lq=1 


© ew ans 
a2 22, 2 BigP pmDantAnn (2.23) 
where ie 
eK, de, SS fl, 
i =o 2 2 Boal 2 B, a3 (p.r, m)L(s, 9, 0)+ 


+-H(q, 8,n)L(p, r, m)—- H(p, r,m)L(q,8,n)— 
a 


—- Hi, 8,n)L(r, p, m)|| (2.24) 
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and the prime denotes that the (m,n)th term of the double summation has 
been taken over to the left-hand side of the equation. Definitions and 


formulae for evaluation of the integrals D 
H(p,r,m) are given in Appendix A. 


pm 


TABLE 1 


Values of the Ay,» coefficients for linear temperature distribution on the border 
strips, and a = 1,o = 0-:73,¢ = 0 








A 5 X10? 108 2°5 X 108 5 x 108 104 
Ai 0°003,747 0°013 55 0°057,08 O'117,3 0°187,5 
Ag, 0°035,20 0°065,34 0°126,9 0°167,0 0°187,4 
Ais —0°000,238 | —0°001,39 | —0°093,73 | —0°004,9 | —0°00I,9 
Ags —0°001,55 | —0°006,45 | —0°008,23 | —0-008,0 | —0°006,3 
As. —0°000,102 0°000,42 0°001,61 0°004,4 O°O10,7 
Aq 0°002,22 0°005,18 0-011,80 0°021,5 0°031,4 
Ai. —9°000,027 | —0°000,10 |—0°000,44 |—0°001,0 | —0-00!,7 
Ass —0°090,24 | —0°'000,49 |—0°001,49 |—0'002,7 | —0°002,6 
As, —0°000,017 | —0°000,05 0°000,07 0°090,2 0°000,4 
Ags —0°000,13 —0°000,31 —0*000,81 —o-001,6 —o°00!,6 
Ass —0°009,011 | —0°000,24 —0°000,35 —0*000,8 —0°001,7 
Ag, 0°000, 34 0°000,67 0°001,83 0°003,9 0°007,1 
Ais —0°000,02 | —0°000,I1 —0°000,4 | —0°000,9 
Ag —0'000,09 | —0°000,33 | —0'000,8 | —o-001,8 
Ase —0°000,03 | —0'000,18 |—0'000,6 | —0'000,8 
Aas —0°000,13 —0°009,26 |—0-009,6 | —0-000,7 
Ass —0°000,06 0°090,04 0°000,1 0°000,3 
Aes 0°900,05 —0°000,22 —0°000,4 —0°000,6 
Ars —0°000,03 —0°000,13 —0°000,5 —0°000,8 
Agi 0°000,13 0°000,35 0°000,8 o-001,8 
Ati —0°000,1 

As» —0°000,4 
Ass —0°000,3 
Ay; —0°000,3 
Ase 0*000,2 
Aes —0°000,I 
Ars —0°000,I 
Ags —0°000,3 
Ags —0°000,3 
Ayo 0°000,6 




















C,,» K(p,r,m), L(p,r,m), and 


The coupled simultaneous algebraic equations (2.21) and (2.23) may be 
solved to any desired degree of accuracy by using an iterative procedure. 
The computational procedure adopted for solving these equations when 
a=1,¢=0,0= 0-73, and A = 5x 10*, 10%, 2-5 108, 5x 108, and 104 
is described in Appendix B. Note that in these solutions it follows from 
symmetry considerations of the basic equations (2.3) and (2.4), that the 
A,,,, can only occur when m-+-n is odd and A,,,, = 0 when m-+n is even. 
Furthermore the B,, ,, can only occur when m-+-n is even. In Tables 1 and 2 
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the calculated A,,,, and B,,,, coefficients are given. It is believed that the 
last figure quoted may be in error by not more than a few units. 


Case II. Border strips being perfect insulators 


This case may be treated in exactly the same way as in Case I. We seek 
solutions of the form 


6 = z s A, q 008 pr€ sin gay 


p=0q=1 


and ib — > > Big X pl€)X{n), 
p=liq=1 


TABLE 2 


Values of the B,,,,, coefficients for linear temperature distribution on the border 
strips, and « = 1, o = 0-73,¢ = 0 








A 5X10? 10° 2°5 X 108 5 x 108 10* 
By 0°267,1 0°519,3 1°159,9 I°919,1 2°885,6 
B's 0°012,6 0°028,6 0°074,4 0°187,4 0°406,7 
Bas 0°001,9 0°008,0 0°031,8 0°059,1 O°051,5 
Bs, O°O12,2 0°027,2 o°061,1 o°118,0 0°200,4 
Bis 00016 0°003,6 0°012,6 0°035,5 O°09I,3 
Bas 0°000,5 o-001,8 0°008,6 0°021,6 0°040,4 
Bs O°001,2 0°002,4 0°006,0 O°O11,9 0°025,5 
By | —0°000,1 | —0'000,4 | —0-002,3 | —0'005,4 | —0-009,6 
Bs 0'001,6 0°003,2 0°008,4 0°018,5 0°035,5 
Bis 0°000,4 0°000,7 0°002,0 0°004,7 0°010,6 
Bes 0°000,1 0°000,2 0°001,6 0°004,6 O’Orl,3 
B35 0*000,2 0°001,0 0°000,3 | —0'000,8 | —0-003,8 
By, 0°000,0 0*000,I 0°000,4 0*001,0 0°003,4 
Bs 0°090,2 0°001,0 O°001,4 0°002,8 0°006,6 
Bes 0°000,0 | —0°000,1 —0°000,3 | —0'000,8 | —0'002,0 
By 0°000,4 0°000,9 O°001,7 0°003,2 0°005,8 
By» 0°003,4 
Bas 0°001,8 
Bs; 0°002,2 
Bas 0°002,2 
Bs —0°000,2 
By, 0°001,8 




















which satisfy the boundary conditions (2.5) and (2.7). Using the appro- 
priate transforms a set of simultaneous algebraic equations defining the 
A,,, and B,,, which are similar to (2.21) and (2.23) may be obtained. Since 
this case is not investigated numerically in the present paper, the actual 
equations will not be given. However, the integrals to be evaluated are 
either those already defined in Appendix A or may be obtained from them 
on integration by parts. 


So far in this section no rigorous mathematical discussion has been given 
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on the validity and convergence of the expansions (2.8) and (2.9). The 
method can only be considered as approximate since the infinite double 
expansions have been taken as finite, although sufficient coefficients have 
been computed in each case to ensure that there is no change in the leading 
terms. Such proofs of validity and convergence are beyond the scope of 
the paper. However, it was thought desirable to check by an independent 
numerical method at least one of the solutions obtained in the present paper. 
Obviously the standard relaxation method is the only available means. 
From Tables 1 and 2 the functions © and % have been evaluated when 
A = 10*, at the nodal points €, 7 = 0-0(0-05)0-3 for +n < 0-3. Hence 
by standard relaxation methods if © and y% are the correct numerical 
approximations to the solutions of (2.3) and (2.4), then the residuals at the 
above nodes must be identically zero, or of such a magnitude that when 
relaxed to zero the changes in the wanted functions are negligible.t (See 
Allen (8).) It was found that the calculated residuals (using mesh length 
equal to 0-05) although being non-zero altered © by less than 0-2 per cent 
and u by less than 0-4 per cent, when relaxed to zero. Spot checks were also 


applied in the region of £, 7 = 0-5 and here the residuals were identically 
zero. 


3. Thermal results 

From Table 1 we can now calculate, for a range of Rayleigh numbers, 
the rate of heat transmission across a cavity having square cross-section, 
vertical hot and cold walls, and a linear temperature distribution on the hori- 
zontal walls or the border strips. As defined in equation (1.8), a dimension- 
less quantity describing the heat transfer is the Nusselt number 


v= (5) 


and on using (2.1), (2.2), and (2.8) we obtain 


fhe) Janine Satta 


p=04 
In Table 3 the calculated values of N are given. These values of V have 
been plotted against G* in Fig. 2. 
Considering first values of A less than 10*, we find from Table 3 that 


N = 145-08 x 10-842, (3.2) 
which is in agreement with Batchelor’s prediction (1.10), where the con- 


+ Note the above statement is correct only if the mesh length is small enough for the 
relaxation solution to converge to a solution of desired accuracy. 
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stant of proportionality was given to be very roughly of the order 10-*. 
For A > 10° the Nusselt number rises rapidly in numerical magnitude until 
in the region A ~ 104, N varies linearly as G? as in Fig. 2. When A =10* 
ordinary conduction across the cavity accounts for 59 per cent of the total 
heat transfer and free convection by laminar flow for the remaining 41 per 
cent. Using the calculated values of N we obtain 


N = 0:16(A/o)* for A ~ 104. (3.3) 


TABLE 3 


Calculated values of the Nusselt number N 








A 5x 107 108 25x 10% 5x 103 104 
Gt 5116 6-084 7°649 9°097 10°818 
N 1'O1 2g I1'O4Ig I*IQ3 I'4I5 1°70, 
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Fic. 2. Variation of heat transfer with Grashof number 





From the limited experimental data of Mull and Reiher (1), Jakob suggests 
that for high aspect ratios N may be represented by the empirical formulae 

N = 0-18(A/o)ta-8® (104 << A < 105). (3.4) 
It is questionable how far beyond the range of the experimental data (3.4) is 
valid, and in particular what is the smallest value of « to be used. Obviously 
not as far down as « = 0 when it would yield an infinite Nusselt number. 
However, (3.3) and (3.4) are in favourable agreement when « = | and it 
seems probable that (3.4) may be valid for a > 1. 

In Figs. 3 and 4 isothermals and streamlines have been drawn for A = 104. 
At this value of A it is apparent that an isothermal core exists in the cavity 
having uniform temperature # = }, also constant vorticity which was veri- 
fied numerically. In the region between the core and the cavity walls, 
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where there exists a continuous boundary layer, @ and y oscillate once near 
the core and then tend smoothly to their appropriate values at the cavity 
walls. The existence of an isothermal core was postulated by Batchelor (3) 
in obtaining the equations describing the flow in the continuous boundary 
layer. 


6=9 0102 03 04 05 0607 0809 10 


S 
| = 


Fic. 3. Isothermals for A = 10* Fic. 4. Streamlines for A = 10* 


























4. Conclusions 
(a) Numerical technique 

In section 2 of this paper a method is developed for obtaining solutions 
of the two-dimensional governing equations in the form of approximate 
double expansions involving orthogonal polynomials. By means of the 
appropriate transforms the governing equations are reduced to the solution 
of two sets of coupled simultaneous algebraic equations. It was found that 
the main computational labour involved was not in solving these simul- 
taneous equations, but in obtaining them. Obviously if only one solution 
of the governing equations was required the labour involved would be a 
serious consideration. In this respect the method given is adaptable to the 
present problem since solutions are required for a range of values of the 
parameters A,¢,anda. It seems likely that the labour involved in obtaining 
solutions by the present method for A > 104anda > 4 would be prohibitive. 


(6) Thermal results 

Preliminary thermal results have been obtained for air enclosed in a 
cavity having square cross-section. The vertical walls are held at different 
temperatures 7, and 7, (7, > 7,), and the temperature along the horizontal 
walls is taken to vary linearly from 7, to 7,. These results indicate that 
for A < 108, N—1 varies as A?, with constant of proportionality 5-08 x 10-°; 
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for A ~ 104, N varies as (A /c)', with constant of proportionality 0-16 (which 
isin favourable agreement with the empirical value 0-18 given by Jakob (2)). 
Finally, when A = 10* an isothermal core exists in the cavity. The vorticity 
in the core has been verified to be constant. In the region between the core 
and the cavity walls there is a continuous boundary layer. 


APPENDIX A 
For completeness we shall give here definitions and formulae for evaluation of 
certain integrals occurring in section 2. For the integrals occurring in (2.21) and 


(2.22) we have 
1 


at, = 2 | X/(€)sinna€ dé = = a {X4(1)(—1)"—X7(0)}, (Al) 
0 
7 1 2m7r _ m_ Yn) 9 
bY, 2/Xx (n) sin mary dn = (may — par r (1M —1)"—X7"(0)}, (A2) 
i 0 


1 
Rir, p,m) X,(£) cos pr sin mr€ dé 


0 


= OniptGm-yp ifm>p 
v r j - 
Ory — Sa ifm < p, (A3) 


and finally , 


T(8,q,n) = [ X,(7n) sin q7rn sin nin dy 


0 
= (q—n)ap_n—(Q+N)ajin, ifg>n 
= (n—q)a,-g—(Q+N)ap+n ifg<n. (A4) 
Again with regard to (2.23) and (2.24) we have: 


1 1 
Dinn = | XmXndE = — [ X;,X de 


0 0 


2) 
=—} > 0b", if m+n is even 
r=1 





= 0, if m+n is odd. (A5) 
1 
Cn ca Xn dg lt ok {X5,(1)— XR (0)}, (A6) 
— 
0 
1 is a] 
K(p,r,m) = [ XI X,X, dé = dj™{X%(1)—X5(0)} — 2 y -  arm™a?, (AT) 
0 n=1 = 
H i xex a as 
i ‘yy ‘a otk... nyu baad ) 
(q,8,n) ea X,X,X, d€ = 2, pat (Xq"(—1)?—X7"(0)} + 5 3 ple tg, 
poi A8 
and — 
1 « ~ 
L(8,q,n) = | XU X,X, dé = > = a"(X"(1)(—1)?— "(0)} + Hy > Pa = ds" a8, 
2 0 tr 
0 p=1 p= st (A9) 


In the above equations (A7)-(A9) the products X,(€)X,(€) and X}(€)X,(€) have been 
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@ 
expanded as cosine series in the range (0,1), obtaining X¥,X, = > d?4 cosnré, 
n=0 


2) 
Xi, X, = > hi cosnré, where 
n=0 
1 
dpa = | XX, df = 8yq 


0 


1 @ 
d?-a — dw? — 2 | X,X,cosnmé dé = 4 > aP(at,_,,—a%,,) forn = 1, 2, 3...., 
r=1 
1 on 
hpa = [ X,X, dé = § > bRat, 
0 n=1 
2 x 
h?a — 2| XX, cosnn€ dé = 3 > b(at_,,—a?%,,) forn = 1, 2, 3.,.... 
a r=1 
The above procedure of evaluation of the required integrals by the aid of Fourier 
series, and in particular the integrals in equations (A7), (A8), and (A9), whilst looking 
rather complicated is in fact much more rapid than the evaluation of formulae 
obtained using the form of X,(€) given by equation (2.11). The Fourier coefficients 
a’, and bf were evaluated for r = 1(1)9 and n = 1(1)60. (Actually af = 0 except 
when r+n is even, bf, = 0 except when r+n is odd). The d?“ and h? are then readily 
obtained by evaluating the above single summations for p,q = 1(1)9,p+q < 10, and 
n = 1(1)25. This task completed, the integrals K(p,r,m), H(p,r,m), and L(p,r,m) 
can then be obtained. The following relationships, obtained on integration by parts, 


K(p,r,m) = —{H(p,r,m)+-H(p,m,r)}, 
L(p,r,m) = —{Lir,p,m)+ L(m, p,r)} 


were found extremely useful in checking the computations. Tables of K(p,r,m), 
H(p,r,m), and L(p,r,m) have been completed for p = 1(1)9, r, m = 1(1)9 for 
r+m < 10. The calculated double array (A5) agreed exactly with that given by 
Young (5). 


APPENDIX B 


The purpose of this appendix is to describe the computational process used in 
solving the coupled simultaneous algebraic equations (2.21) and (2.23) when a = I, 
¢é = 0 for various Rayleigh numbers. This process could be used, without modifica- 
tion, to treat other values of the parameters a and ¢. 

Some mention should first be made of the way in which an individual equation, 
say for Aj», is recorded. The first double summation on the right-hand side of (2.21) 
will involve, by virtue of the symmetry properties of the functions X,(£), the B,,, for 
r, 8 = 2(2)8,r+s < 10. A double array of the product a§ bj is then calculated from 
values already computed using equations (Al) and (A2). The quadruple summation 
Yi, is treated in a similar fashion. First when r and s are both even, taking r,s = 2(2)8, 
r+s < 10, double arrays for 


pR(r,p, 1) T(s, 9g, 2) —qR(s, q, 2) T(r, p, 1) (Bl) 


were evaluated using (A3) and (A4) for p = 1(2)9, q = 2(2)8, and p+q < 10. 
Secondly, when r and s are both odd (i.e. r, 8 = 1(2)9, r+8 < 10) double arrays of 


(B1) were evaluated for p = 2(2)8 and g = 1(2)9. Thus given the values of the A, , 
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for A ~ 104, N varies as (A/c)*, with constant of proportionality 0-16 (which 
isin favourable agreement with the empirical value 0-18 given by Jakob (2)). 
Finally, when A = 10‘ an isothermal core exists in the cavity. The vorticity 
in the core has been verified to be constant. In the region between the core 
and the cavity walls there is a continuous boundary layer. 


APPENDIX A 
For completeness we shall give here definitions and formulae for evaluation of 
certain integrals occurring in section 2. For the integrals occurring in (2.21) and 


(2.22) we have 
1 


2nar 
as a X,(&) sin nwé dé = . aes (1)(—1)"—X2(0)}, (Al) 
r 1 2m7r oa . aie ‘ 
= 2 Z| X}(y) sin my dn = Gomer yt ps eC -1)™—X""(0)}, (A2) 
1 
R(r,p,m) ( X,(€) cos pr€ sin mre dé 
d 
an +p Tv oe if m > Pp 
On+p—2p—m ifm < P> (A3) 


and finally 1 


T(s,q,n) = X,(7) sin q7rn sin nn dy 
0 


= (q—n)ap_n—(q+n)ajin ifg>n 


= (n—q)a_g—(q+n)ajen, ifg <n. (A4) 
Again with regard to (2.23) and (2.24) we have: 
1 1 
Dn = { XXndé = — [ Xj,Xp dé 
fy) 0 
4 
= —})> 0b", if m+n is even 
r=1 
= 0, if m+n is odd. (A5) 
, 
C, = X,, dé = —{X™(1)—X”(0)}, (A6) 
0 
1 ies saree 1 
K(p, r,m) = {xy X,Xm dé war am X5(1)—X5(0)} He ed a dy; mans (A7) 
0 *, n=1 ie 
t 2 
H ,8, “a rey jie Fim the cy Ha " 
(q,8,”) ) X, X,X,, dé 2 mm hg x”"(—1)?—X7"(0)} 4-8 >> mal ug, 
and $3 a“ 
1 ao cs 
L(8,q,n) = | X,X,X, dé = > 4.2 ds"{X"(1)(—1)?— X70 Ha x @. 
} fey n é < pin? ? t gf )( ) XG t+ a) ~ p saps da” ay F 


p=1 
(A9) 
In the above equations (A7)-(A9) the products X,(€)X,(£) and X}(€)X,(€) have been 








HEAT TRANSFER BY LAMINAR FREE CONVECTION 271 


@ 
expanded as cosine series in the range (0,1), obtaining X,X, = > d?“cosnzé, 
n=0 


@ 
Xi, X, = > kh*cosnn€, where 
n=0 


a 
dg¢ = | X,X, df = 5, 


0 


1 
d?-a — dP = 2 f X,X,cosnmé dé = } s a?(a?,_»;—a%,,) forn = 1, 2, 3.,..., 
r=1 
1 : a 
hpa = [ X,X, dé = 4 > bag, 
0 n=1 
1 2) 
hea — 2 -- X,cosnmé d& = 4 > bP(at,_,,—a?,,) forn = 1, 2, 3,.... 
0 r=1 
The above procedure of evaluation of the required integrals by the aid of Fourier 
series, and in particular the integrals in equations (A7), (A8), and (A9), whilst looking 
rather complicated is in fact much more rapid than the evaluation of formulae 
obtained using the form of X,(€) given by equation (2.11). The Fourier coefficients 
a’, and bf were evaluated for r = 1(1)9 and n = 1(1)60. (Actually aj, = 0 except 
when r+-n is even, bf, = 0 except when r+-n is odd). The d?¢ and h?- are then readily 
obtained by evaluating the above single summations for p,q = 1(1)9,p+q < 10, and 
n = 1(1)25. This task completed, the integrals K(p,r,m), H(p,r,m), and L(p,r,m) 
can then be obtained. The following relationships, obtained on integration by parts, 


K(p,r,m) = —{H(p,r,m)+H(p,m,r)}, 
L(p,r,m) = —{L(r,p,m)+ L(m, p,r)} 


were found extremely useful in checking the computations. Tables of K(p,r,m), 
H(p,r,m), and L(p,r,m) have been completed for p = 1(1)9, r, m = 1(1)9 for 
r+m < 10. The calculated double array (A5) agreed exactly with that given by 
Young (5). 


APPENDIX B 


The purpose of this appendix is to describe the computational process used in 
solving the coupled simultaneous algebraic equations (2.21) and (2.23) when a = I, 
¢ = 0 for various Rayleigh numbers. This process could be used, without modifica- 
tion, to treat other values of the parameters a and ¢. 

Some mention should first be made of the way in which an individual equation, 
say for A», is recorded. The first double summation on the right-hand side of (2.21) 
will involve, by virtue of the symmetry properties of the functions X,(£), the B,, for 
r, 8 2(2)8, r+s < 10. A double array of the product a§b{ is then calculated from 
values already computed using equations (Al) and (A2). The quadruple summation 
Yi. is treated in a similar fashion. First when r and s are both even, taking r,s = 2(2)8, 
r+s < 10, double arrays for 


pRir, p, 1) T(s,q, 2)—qR(s, q, 2) T(r, p, 1) (Bl) 
were evaluated using (A3) and (A4) for p = 1(2)9, gq = 2(2)8, and p+q < 10. 
Secondly, when r and s are both odd (i.e. r, 8 = 1(2)9, r+8 < 10) double arrays of 
(B1) were evaluated for p = 2(2)8 and g = 1(2)9. Thus given the values of the Ay, , 
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and B,,,, anew estimate of A, ,may be obtained using equat ion (2.21). The first doub| 
summation involving the B,,, coefficients may be readily obtained. The ee le 


summation y; is written as 
“E22 ail. 23P.9)> (B2) 


p=1 g=1 


y(1,2;p,q) = s s B, {pRr,p, 1)7(s,q, 2)—qR(s,q, 2)T(r, p, 1)}. 
r=l1 8=1 


where 


From the double arrays (B1) calculated for p = 1(2)9, g = 2(2)8, and p+q < 10 
we can then calculate the required y(1, 2; p,q) by performing double summations over 
the B,,. Finally, by means of a double summation over the A,,, and the y(1, 2; p,q) 
as in (B2), y,,. is obtained. 

The iterative procedure adopted for the solution of (2.21) and (2.23) was as follows. 
Initial estimates of the B,,,, for m,n = 1(2)9, m+n < 10, were obtained using the 
equation 


1 2 re 
(- =m +3 5 Dam Dant tt) Bu. = i AC, Ca— 5" BygDomD 


p=1 @=1 


the quadruple summations y,,, being neglected at this stage as they are small com- 
pared with the term involving the double summation in A,,,, about which at this 
stage we have no information. Returning to equation (2.21) we can now calculate 
the A,,,, for m even and n odd, using 


1 1 ® @ 
(4 mttnt)A mn = 5 :> > B, a}, jn) 


r=1s=1 
SINCE Ym,» = 0. Values of 


> > (4 p,q cos ar br —isingay on) 


p=1 q=1 


can now be calculated for m and n even and so we obtain the first estimates of B,,, ,, 
for mand n even. A fairly accurate set of B,, , now exist and from these the y,,,, may 
be obtained giving a complete set of the A, 4; Am,, is then computed and then the 
current set of B,,,. Fortunately after a few complete iterations, the values of the 
quadruple summations y,,, and A,,, remain steady and the main labour involved 
in each iterative step lies in controlling the sensitive relationship between the main 
By, and A», » coefficients. The higher the value of A the more serious is the situation. 
This difficulty was resolved by setting up new equations for the A, , (takingm+n = 3 
and 5) and B,,, (taking m+n = 2, 4, and 6) which involved only the coefficients 
indicated, the remaining coefficients having been added in numerically as these 
achieve their final values at an early stage in the iteration. It was then much easier 
to control those few coefficients which oscillated about their true values. However, 
even in the case A = 10‘, no more than five complete iterative steps were required. 
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SUMMARY 
The case of flow at large Reynolds number is investigated, special attention being 
paid tc the problem of finding deviations from the well-known inviscid flow. The 


fundamental equations in this case are reduced to a single second-order differential 
equation of the singular perturbation type, which is quite similar to the corresponding 
equation obtained in the case of two-dimensional radial flow. Utilizing this similarity 
and the knowledge obtained from the solution of the two-dimensional case as a guide, 
the solutions which contain shock waves are studied in detail. 


1. Introduction 


THE purpose of the present paper is to investigate the solution of the Navier— 
Stokes equations in the case of steady, three-dimensional, radial flow of both 
source and sink types, taking into account the effects of compressibility, 
viscosity, and heat conduction. This problem in the absence of viscosity 
and heat conduction (inviscid problem) is one of the few examples whose 
exact solutions are known and expressed in simple closed form. Thus it is 
of special interest to find the deviation from the inviscid solution due to the 
effects of viscosity and heat conduction. 

The same kind of problem has been discussed in the corresponding case 
of two-dimensional radial flow. The present author (1) reduced the problem 
to a differential equation of the first order for source flow and solved it for 
a certain Reynolds number. Levey (2) investigated the equation in detail 
and proved the existence of shock waves in a class of flows at large Reynolds 
numbers. The case of sink type flow was considered by Wu (3). He ex- 
pressed a solution in an analytical form at large Reynolds numbers using 
Lighthill’s method. Prosnak (4) investigated the problem for both cases of 
source and sink flows from a different point of view in order to find the 
separate effect of either viscosity or heat conduction alone. 

Under the assumption of large Reynolds number, the system of equations 
for the present problem is reduced and simplified to one equation: 

d*w 


dw 2 
_aatbese — 2 — = >] === s] —_— 


+ Now at Tokyo Electrical Engineering College, Kanda, Tokyo. 
(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 3, 1958] 
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where w = u/e,,x =17,/r, B = (y—1)/(y+1) (see (8), (9), (15), (25) below) 
and a is a small number proportional to the reciprocal of the Reynolds 
number. The corresponding equation for the two-dimensional case is found 
in the papers referred to above, and is 
/ 
ay) ’ 

aut — (1—ws) + w(1—Bw?) = 0, t = logr. (B) 
Comparing these two equations, we can see that they are similar except 
for the factor 2/x in the equation (A) and solutions of these equations 
would be expected to be similar. Actually, the factor 2/z is almost constant 
in a narrow range of 2, and thus a sclution of (B) which includes shock 
waves may also be expected to be a solution of (A), since in the shock wave 
the change in w takes place in a narrow region of x. 

The equation (A) is non-linear and of singular perturbation type, having 
a small coefficient a of the highest derivative. It is difficult to express the 
solution in a single analytical form in the whole range of x. Utilizing the 
similarity between (A) and (B), and using the knowledge of the solution 
of (B) as a guide, the domain of x is divided into several sub-regions where 
approximate solutions can be obtained. By connecting these approximate 
solutions in different regions using boundary layer technique (Carrier, 5), 
the solution can be constructed. 

After reducing the fundamental equations to the equation (A) in section 2, 
the properties of the solution of (B) are displayed in section 3 for the purpose 
of obtaining qualitative information avout the solution of (A). Detailed 
study is given in section 4 to the typical solution of (A) which, in the 
limit of large Reynolds number, reduces to a solution consisting of inviscid 
solutions and shock waves. Other types of solutions are briefly discussed 
in section 5. 


2. The fundamental equations 
Let r, u, p, p, T denote respectively the radial distance from the origin, 

velocity, pressure, density, and absolute temperature. The system of 
Navier-Stokes equations specialized to the case of three-dimensional steady 
radial flow of a perfect gas with viscosity 4 and heat conductivity A are as 
follows: 
The continuity equation is pur? = m, (1) 
where m is a constant and 

m > 0 for asource flow, 

m <0 for asink flow; 
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the momentum equation is 
du _ 4p hk 2,)\ oS. 
pu = 5 5 (ru)| —— 
dr adr'3 dr\" dr’ . 
the energy equation, after one integration and ae the integration con- 
stant as m£, is 


dT 
2me,, T+-(m+§uT )u? = Dr? + Sur? 


(2) 


du? 
dr 


these equations are supplemented by the equation of state 
p = pRT. (4) 
The equations for the inviscid fluid are obtained by putting A, » = 0 in 
the above system of equations. Then (2), (3) become 


+2mE; (3) 


du dp 
sees SEE: inmates 5 
ew ar dr’ “) 
2c, T+u* = 2E = 2,7, (6) 


where 7), is the temperature at the stagnation point. Eliminating p, p, 7 
from (1), (4), (5), (6), we get an equation for wu, 

duly 2 r 
ee §— 2, T,) += u(u*—2e, T,) = 0, (7) 


where y is the ratio of the specific heats. Introducing a new variable w 
defined by 


w= u/c, (8) 

where c, is the critical velocity defined by 
c} = 2Bc,, 3 B — (y—1)/(y+)), (9) 
(7) is reduced to a (1 —w) +? (1 — But) = 0. (10) 


This equation gives, on integration, 


Fa wri( But) oO 


; (11) 


where we have put r = r, at w = 1. The graph of w/(r) given by (11) is 
plotted in Fig. 1. The graph shows almost the same features as in the two- 
dimensional case, i.e. it has two branches—one supersonic one subsonic— 
inr > r,, terminating at r = r, with an infinite gradient and there is no 
solution in r < r,. 

For r > r, the slope of the curve w(r) of (11) is much smaller than unity, 
and consequently viscous effects should become comparatively unimpor- 
tant there, and so the inviscid solution (11) itself can be considered to 
be an approximate solution of the system of equations (1)-(4) with non- 
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zero but small A, », provided E = 2c, T, and r is large enough. Since our 
purpose is to find the deviation from the inviscid solution (11) due to the 
effects of heat conduction and viscosity, it is natural to seek a solution 
which approximates to the inviscid solution (11) for large value of r. More 
precisely, we shall find a solution which reduces to the inviscid solution 
(11) when A, « + 0 at large but finite r. Accordingly EZ is determined as 
E = 2c,,T). We must distinguish this 
solution from a similar but different + 
one which reduces to the inviscid 
solution as r > oo when A, » are small 

but finite. The second solution cannot 1 Pi. 
always be found, because the super- 
sonic branch of (11) does not satisfy 
(3) when r > oo if A, pw are finite and 
fixed. 

In the papers quoted above, the effects of varying A, » were estimated 
to be inconsiderable in the case of two-dimensional flow. For the purpose 
of making the problem as simple as possible, we shall confine our attention 
to the case when A, p» are constants. 

Now, eliminating p and p from equations (1)—(4) and introducing a non- 
dimensional variable z, 











' 
' 
' 
‘ 
afl, 


0 A »T 
Fia. 1 





r _ 3m {z>0 for a source flow 


= — © at a 2 
oes 4u \|z<0 fora sink flow, (12) 
the equations for u, 7’, assuming A, py, c, constant, become 
d*u du y—l{d(c,T\ 2c,T 
Ge 4:tOy-- Pi 96-00 Bo an 1 
dt | : a . y Fal u ) zu (13) 
A ee 2 d (u? 
Cy T+(1+22z)> m 8 qa? T)+z i) tT (14) 


where o is the Prandtl number and is also assumed to be constant. Unlike 
the case of two-dimensional flow, r* in (12) is not a non-dimensional con- 
stant interpreted as a Reynolds number but has the dimension of length. 
Since equations (13) and (14) do not depend explicitly on r*, the condition 
for (1)-(4) when A, x — 0 at finite and fixed r becomes the condition at z > 0 
for (13), (14). It is convenient to introduce a new variable y and a small 
number a such that 

pe as {x >0 for a source flow 
“ : r* \a<0 fora sink flow 
and a is interpreted as the reciprocal of the Reynolds number Re (a-! = Re) 
and assumed to be small. 


(15) 
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Using (15) together with the non-dimensional variables 


T : 
nem. Pe (16) 
Cy 7, 
where 7; is the value of 7 at r = r, in the inviscid solution, equations 
(13), (14) are transformed into 


dw. dv 4. a {0\ 28) : 

— —_- — WA, = = a Daa as : l 
+r dy? i+ (2a ad iz a yldy\w/ yu (17) 
ste tt e+ Is 


The problem is thus reduced to finding a solution of (17), (18) which reduces 
to the inviscid solution when «> 0 at finite and fixed y. By a further 
transformation x = ly =n,/r, (19) 


equations (17), (18) become 


dw dw 2%  £élfd/@\ 26) 9 
Gat ge a’ = — ala) +3 ap (20) 
= 2a) yt _ 3d) y—1 dw"\ 
+ l SP ' 2 
e+ | ss =|- 2 of at 2 dx (21) 


which are similar to the corresponding equations in the two-dimensional 
case. 

Now let us write (21) in the different form 
4o a 1 tye tt 


045 


a) 
where A is an arbitrary constant. The last term is small if 1—(40/3) (= e, 
say) is small enough, i.e. if o ~ 0-75 (for air o = 0-72). Now 2 exceeds a 
certain fixed value x, (say) (we are considering the solution in a range 
r < R fixed) and it is likely that w will not exceed the maximum value B-+ 
in the inviscid solution. Then we have 


yond A aan f ctanel (149) 4 Nu dx < ae “(e+ =)B. 1 


—14 1 2 
-_ Ae terre _ 1 pectates [tome (1 — 2) 4 Sls dx, (22) 





es 


With regard to the term Ae~“4*3~", 4 must be zero for « < 0, otherwise it 
becomes infinite as a > 0 for finite 2; on the other hand, for « > 0, Ae-4030 
becomes smaller than any finite power of «. Then (22) becomes (for both 


40 y—1 1 2 
i > wt = ofe+>). (23) 


cases a 2 0) 


3 
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Substituting (23) into (20) and neglecting small order terms, we get finally 


d*w dw 2 
on ef) eg, cath — teh oe 24 
ut (ws) 52 + = wo(1— Bt) = 0, (24) 
». 
where —Y «=a. (25) 
yt+l 


If we put a = 0 in (24), the equation reduces to equation (10) having the 
solution (11). 


3. Qualitative properties of the flow from the solution of the two- 
dimensional case 


Let us compare equation (24) with the corresponding equation in the 
two-dimensional case, namely, 
d*w dw 


aw* , (1—w!) + (1—Bu*)w = 0, t = logr, (26) 


neglecting small quantities. Since (24) and (26) are similar, except for the 
extra factor 2/2 in the last term of (24), it is to be expected that the solutions 


ad 


Ww 














0 nr 7 7 
Fic. 2. Sketch of a solution curve for the two-dimensional 
source type flow (a > 0) 


will be similar, the independent variables being x in (24) and ¢ in (26). In 
particular, the local properties of the solution of (24) will be the same as 
for (26), since the factor 2/x in (24) veries little if x ranges over a small 
region. We shall use this fact to find the qualitative properties of the 
solution of (24). Before doing this, we shall consider the general properties 
of the solution of (26). 

The solution of (26) has not been found in finite form since it is of singular 
perturbation type. But the important characteristics of the solution are 
known and may be found in the papers mentioned above. The typical 
solution curve for a > 0 is shown in Fig. 2. 

The curve starting from the stagnation point at r = r,, (1) rises and 
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(2) passes the sonic line w = 1, then (3) it goes along the supersonic branch 
of the inviscid solution w, until near r = ro, (4) it falls down to the subsonic 
branch of another inviscid solution w, and finally (5) approaches the 
stagnation point at infinity along wy. At the limit of a > 0, parts (3) and 
(5) approach the corresponding parts of inviscid curves while (4) approaches 
the shock wave. The parts (1), (2) approach the line P’r, (‘negative shock’). 
Since the flow is not isentropic but is affected by viscosity and heat con- 
duction, it is natural that the solution as a > 0 corresponds to different 
inviscid solutions which have different entropies. There is an infinite 
number of solutions of this type depending on the different choice of 




















pi pie 
“i, = 
Pd ‘ /“o 
Pf2)4 P'! 
his ail a Aes 
. t 
1 iN i\ 
PN 9 ae 
— ih, at a 
(1) (ii) 


Fic. 3 


11,17, %,- Two of them are, however, independent. For example, if r,, 7, are 
given, then r, must be determined so as to satisfy the shock condition at rg. 

Now, if we move r, towards r, which is fixed, r,; moves also towards r,, 
and part (3) of the curves becomes smaller and finally disappears. 

Special consideration must be given to this case, or more precisely, to 
the case where (r,/r,;)—1 < O(a'), and the curve illustrating the solution 
is as shown in Fig. 3 (i). 

The solution in the transition regions (2), (4) is described by a kind of 
transonic equation. But such a region is small (O(a!)) and becomes zero 
when a0 and then the solution at the limit consists of two parts: 
(1) negative shock, and (5) the subsonic branch of the inviscid solution wy. 
It must be remarked that the limit of this solution corresponds to a single 
inviscid solution. 

Solutions of the same type which consist of the two parts (1), (5) shown 
in Fig. 3 (ii), are possible, though the transition from (1) to (5) is different 
from that in the solution shown in Fig. 3 (i). 

These main properties of the solution for a > 0 appear also in the case 
of a < 0 (sink type flow) except for the fact that the roles of the subsonic 
and supersonic branches of the inviscid solution are interchanged. A typical 
solution corresponding to that for a > 0 in Fig. 2 is illustrated in Fig. 4. 

The curve starts from maximum velocity B- as r > 00 on the supersonic 
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branch of the inviscid solution w, goes down to (3) (the subsonic branch of 
w ) through (4) and thereafter passes parts (2), (1) and finally ends at the 
maximum velocity w = p>. 

In the limit a + 0, the parts (5), (3) approach to the inviscid solutions 
and the part (4) approaches the shock wave. Since shock transition takes 
place from low to higher values of the entropy and w, has higher entropy 
than wo, the transition must be from w, to wy as was true for the case 
a > 0. Solutions of a special type similar to that shown for the case a > 0 
in Fig. 3 also exist. 


Bp? 














0 








Fic. 4. Sketch of a solution curve for the two-dimensional 
sink type flow (a < 0) 


It is worth while noting that equation (26) resembles in many ways 

van der Pol’s equation for relaxation oscillations, viz. 
d*y 
dt? 
For example, the rapid transition occurring between two regions (2) and 
(3) or (3) and (4) in Figs. 2 and 4 is due to the same causes which induce 
the well-known jerky wave form in relaxation oscillations. 

Returning to the present problem of finding the solution of (24) there 
might be a solution which includes transition regions something like (1), (2), 
(4) since the regions are so narrow that the factor 2/x in (24) has little 
influence on the solution there, while the solution of (26) near the regions 
(3), (5) is essentially net so different from the inviscid solution; the same 
situation is also to be expected in the solution of (24). These points will 
be clarified in the following sections. 


d 
agi (ly) +y = 0. 


4. Typical solution of equation (24) 

Corresponding to the typical solution of (26) which was described in the 
last section, there is a similar solution of (24). This is derived by sub- 
dividing the range of x and finding approximate solutions in the sub-regions. 
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(2) passes the sonic line w = 1, then (3) it goes along the supersonic branch 
of the inviscid solution w until near r = ro, (4) it falls down to the subsonic 
branch of another inviscid solution wy, and finally (5) approaches the 
stagnation point at infinity along w). At the limit of a > 0, parts (3) and 
(5) approach the corresponding parts of inviscid curves while (4) approaches 
the shock wave. The parts (1), (2) approach the line P’r, (‘negative shock’). 
Since the flow is not isentropic but is affected by viscosity and heat con- 
duction, it is natural that the solution as a > 0 corresponds to different 
inviscid solutions which have different entropies. There is an infinite 
number of solutions of this type depending on the different choice of 
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11,7, %- Two of them are, however, independent. For example, if r,, 7, are 
given, then r, must be determined so as to satisfy the shock condition at ro. 

Now, if we move r, towards r, which is fixed, r,; moves also towards r,, 
and part (3) of the curves becomes smaller and finally disappears. 

Special consideration must be given te this case, or more precisely, to 
the case where (r,/r,)—1 < O(a'), and the curve illustrating the solution 
is as shown in Fig. 3 (i). 

The solution in the transition regions (2), (4) is described by a kind of 
transonic equation. But such a region is small (O(a?)) and becomes zero 
when a0 and then the solution at the limit consists of two parts: 
(1) negative shock, and (5) the subsonic branch of the inviscid solution wy. 
It must be remarked that the limit of this solution corresponds to a single 
inviscid solution. 

Solutions of the same type which consist of the two parts (1), (5) shown 
in Fig. 3 (ii), are possible, though the transition from (1) to (5) is different 
from that in the solution shown in Fig. 3 (i). 

These main properties of the solution for a > 0 appear also in the case 
of a < 0 (sink type flow) except for the fact that the roles of the subsonic 
and supersonic branches of the inviscid solution are interchanged. A typical 
solution corresponding to that for a > 0 in Fig. 2 is illustrated in Fig. 4. 

The curve starts from maximum velocity B-+ as r -> co on the supersonic 
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branch of the inviscid solution w4, goes down to (3) (the subsonic branch of 
w») through (4) and thereafter passes parts (2), (1) and finally ends at the 
maximum velocity w = B-+. 

In the limit a > 0, the parts (5), (3) approach to the inviscid solutions 
and the part (4) approaches the shock wave. Since shock transition takes 
place from low to higher values of the entropy and w, has higher entropy 
than wo, the transition must be from wy to w, as was true for the case 


a > 0, Solutions of a special type similar to that shown for the case a > 0 
in Fig. 3 also exist. 


Bp? 


Ww 














0 








Fic. 4. Sketch of a solution curve for the two-dimensional 
sink type flow (a < 0) 


It is worth while noting that equation (26) resembles in many ways 
van der Pol’s equation for relaxation oscillations, viz. 

d*y 
de 
For example, the rapid transition occurring between two regions (2) and 
(3) or (3) and (4) in Figs. 2 and 4 is due to the same causes which induce 
the well-known jerky wave form in relaxation oscillations. 

Returning to the present problem of finding the solution of (24) there 
might be a solution which includes transition regions something like (1), (2), 
(4) since the regions are so narrow that the factor 2/z in (24) has little 
influence on the solution there, while the solution of (26) near the regions 
(3), (5) is essentially net so different from the inviscid solution; the same 
situation is also to be expected in the solution of (24). These points will 
be clarified in the following sections. 


d 
agi — (vty = 0. 


4. Typical solution of equation (24) 

Corresponding to the typical solution of (26) which was described in the 
last section, there is a similar solution of (24). This is derived by sub- 
dividing the range of z and finding approximate solutions in the sub-regions. 
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Detailed consideration is given mainly to the case of source type flow 
(a > 0) since the second case (a < 0) is similar to the first case, at least 
qualitatively. 

In the regions (3), (5), the solutions of (26) are not so different from the 
inviscid solution and the term aw?(d?w/dt®) can be considered to be small 
there, since for the inviscid solution d*w/dt? remains of the same order as 
the other terms except in the neighbourhood of t = 0. In equation (24), too, 
the inviscid solution wo, satisfying the equation 


2\dw, , 2 2 97 
— (13) 50 += wy(—Bres) = 0, (27) 
— Rap2\(1—By/48 
and given by a a ual? ee) , (28) 


is considered to be a good approximation for w in (24) when a is small so 
long as aw*(d*w/dx?) remains small. On the other hand, in the regions 
(1), (2), (4) the term aw?(d?w/dx?) cannot be neglected and wy ceases to be 
an approximation to w; different approximations are then needed. Fitting 
together these approximate solutions in the different regions, we finally 
obtain the overall solution. 

We begin by finding the deviation of w from wp», necessary to fit the 
approximate solution in (3), (5) (from now on, for convenience, we shall 
use these notations also for the corresponding regions in the solution of 
(24)) to the approximate solutions in the neighbouring regions (1), (2), (4). 
For this purpose, we put w= wtf (29) 


where we suppose Sf <u. 
Substituting from (29) to (24) and neglecting small quantities of order 
higher than f?, we get 
mee fs we, lw, 1 1—3Bw? 
mote k Uf 4. Sed aeulieneete ae ae ane : 
f as \s ii (+2 W ax wi \s _ (30) 


where primes denote differentiation with respect to z. By a transformation 


f — e*'%q, 
r 
where Zz =e. . 1) dx 31 
me 5 Bee ae eh) 
0 
equation (30) is transformed into the standard form 
9" —a-*y(x)g = —whe-#, (32) 
where x(x) = xo(x)+ax,(x)+a?x,(2), 


9 Ps 


1/1 8 % 4 % 
ue) =3{a-)>  e)= (4), yey = 28%. a 
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Since equation (32) is of Green’s type, its asymptotic solution for small 
values of a can be obtained by Jeffreys’s method (6). In the region where 
Xo> Xi X2 Temain of comparable orders of magnitude, two independent 
solutions g,, J, of the homogeneous part of (32) are, to a first approximation, 

iven by ne 
g . 91, J2™~ x *e*8, 


zx 
where 7 = | x? dx. (34) 


Using the fact that the Wronskian g; g,—g, 92 is equal to 2/a, the approxi- 
mate expression for f becomes 


f — et/ag 
~~ x7 | Ac#+ma+ Bele-wia_ta e@+na ( wy x teeta dx + 
+a e@-n/a ( ws x te-e-me dz}, (35) 
where A, B are arbitrary constants. 
We consider firstly the case a > 0 (source type flow). Since w, consists 


of two branches, wy < 1 (subsonic branch) and w, > 1 (supersonic branch), 
we must consider the solution for each branch separately. For w, < 1 we 


have me r 
1/1 aw, ({, 4Bw 
n= | x! dx = | ea ) —_ er . | | dx 


= z—alog(1—w?)-4(1—Bu?)+..., 
Se han a ee 


* es a 2 
ws Wy «xi\ur 








and (35) reduces to 
f ~ wo| V2 A(1—Bo§)4e% + V2 B(L—w§)-( 1 — Bes) — 
—a(1—oo§)4e [05 wo(1—we5)-1(1 ge dar + 
+-a(1—w3)-"(1—Bweg) {we wo(1—Bu§) de +...}. (37) 
The third term can be shown to be small by integrating by parts: 
a(1—Bueg)te%* fre w9(1—w§)-4(1— Brose dx 
= —a*wp w3(1—w?)-*+-a?(1—Bw§)-e** x 
x [etn ugh —w)-401—Bud)} de. (38) 
The behaviours of wy and z for small x are 


1 


1 
W, ~ ke’, tw — sa 


where k = (1—p)4-P28, 
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so that behaviour of f/w for small x becomes, from (37), 

flwy ~ V2 Ae~VGak*2) 4.4/2 B+ fak*2', 
Now f must tend to zero as a > 0 at a fixed x (say 2), and so B+ 0 as 
a0. Finally, we get 
fiwy ~ V2 A(1—Bw?)-te*'4+-a(1 —w§) (1 — Bug) | Wo Wo(1—Bw§)-? dx + 

4-2 B(1—w?)-(1—Bw?)+.... (39) 

In the first term e2** is of a higher order of small quantities than any finite 
power of a and thus the first term will be quite small. Accordingly, for 
the limiting case when a is small, the deviation, f, of w from wy is also small 
(f = O(a)) in the region (5) (a > 0, w) < 1). But it is possible that the 
term Ae** may become large at a certain point. Introducing z, instead 
of A, where z, is given by V2A = e-2/a, (40) 


the first term becomes (1—Bw?)-1e%—2a, 


Now, if it happens that z = z, at = 2, (say), the magnitude of this term 
will vary suddenly from a small to a very large value, as x passes through 2» 
(r decreasing), since a is assumed to be a very small number. This feature 
may be incorporated in the following approximate formula valid near 
x = 2, for (39). From (31) we have 


2(z—2) = | (7 y dz ~ (ie 1) @—a9), 
0 01 


assuming that x—z, is small and denoting by w,, the value of wy, at 2p. 
Thus, (39) becomes 


f ~ toy(1— Buch) Aelndd-a¥e-20!6 (w < a), (41) 
neglecting small quantities. 
If we change the variable x to according to the relation 


x = 4 +a€, (42) 
(41) gives - f ~ wo, (1—Bw3,)“relmed-nE (€ < 0). (43) 
The choice of z) or x) (which corresponds to r, in Fig. 2) is completely 
arbitrary so long as Xo(%»), x(%»), X2(%p) in (33) are of order unity. This 
means that such a sudden change can happen anywhere in the region 
1 > x > 2, excluding x, ~ 1. 

Having now determined a position of sudden change x = x», we proceed 
to find an approximate solution which is valid in a narrow region (4) near 
« = 2. Changing the variable from z to é by (42), equation (24) becomes 

lw 
wit z = i w(1—Bw*) = 0. 


eat a tad 
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Neglecting the last term since it is multiplied by a and is small provided 
aé remains small, we have the approximate equation 


d?w dw 
i... to o- 
won ( wy 0, (44) 
which on integration gives 
dw 1 
Es” = ¢, (45) 


where c is a constant of integration. Choosing c so that (45) assumes the 
form 


dw 1 
dé = =, (wW— wy) wae ); (46) 
where wW,+W, = C¢, Uyw,=1 (uw, <1), (47) 


and integrating (46) in the region w, > w > w,, we have 
Ww, w , 
——1_ log(w—w,) — ——2—_ log(w,—w) = €+¢’, 48 
w,—W, g( 1) w,—W, g(w,—w) = §+ (48) 


where c’ is an integration constant. The integration constants c, c’ (or w,, c’) 
can be determined by comparing equation (48) at large negative ¢ (where 
aé remains of order unity) with equation (43). 

Considering the behaviour of (48) as £ - —oo, we have 


w~ w,+(w,—w,)"* Wig Wa—wi)e’/wr p(wa—wi kw 








1—wu? Vut 2 2 ¢)2 
= WwW (1—wi)e’/wi p(1—wpé/wi 
+( w, ) e € ‘ (49) 
Comparing (43) and (49), we obtain the values of w, and c’, viz. 
UW, = Vow (50) 
1—wuj 1w} 2) 47 Han? 
Wy EA UNE NY == Woy(1— Bw, )-*. (51) 


It is worth noting that (46), (47), (48) are the same as the equations for the 
velocity distribution deduced from the Navier-Stokes equation for the case 
of plane shock waves and that the second condition in (47) corresponds to 
Prandtl’s relation. Thus, the transition near x = 2, is seen to have the 
same feature as that in plane shock waves. In the case of two-dimensional 
source type flow, the presence of such a shock wave was derived by Levey (2) 
using quite a different method. 

Now the solution (48) must be fitted to a supersonic branch of another 
inviscid solution for large positive values of €. At large £ > 0, (48) becomes 





2 \ wk, 
ne aa a yr ett wig, (52) 
Wy = Wg,“ \1 wh, 
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where we have used (50), (51). On the supersonic branch (wy > 1) equation 
(35) for f reduces to a slightly different form of (39) and 
f/wy ~ V2 A'(1—Bw?)(w3 —1)-1+- v2 B’(1—Busg)*e* 
Fall Bink Ba ns 
tect —Bw%) | wo wo(1—Bws) da, (53) 
with arbitrary constants A’, B’. Putting 
\2 B’ = B"e-20a (54) 
and denoting by wp: the value of wy (supersonic branch) at x = x, we get 
an expression for f valid near x = 2% (2 > 2p): 
‘ ade ; ee Se. sie 
f ~ V2 A’ wo2(1—Bu2.)(w2.—1)-!+ B’w(1—Bwg.)te-“A ree", (55) 
neglecting small quantities. Comparing the expression (55) with (52) we get 
l , 1—uj, \wht lie 
Wo = We = —, BR’ = a1 , A’=o0/(l1). (56) 
Wo1 —Bui 01 
Since the term B’e* becomes quite small when O(2—x,) becomes larger 
than a, the solution w in the region (3) is again approximated by the 
inviscid solution wy, but this function w, does not correspond to the same 
inviscid solution used before (which was of the form given by (28)), but 
has a different form given by 
1 — Bw2\t-By48 r, st 
vk = wil? B ° k= +, (57) 
1-8 ry 
corresponding to a different entropy. Since at x = 2», wy in (28) is wy, (< 1) 
and wy in (57) is wo, (> 1), we can determine & by using the first condition 


of (56), whence 7 % (i —B/uy, i Byap (58) 
ry Woy 1—Buyy 


It should be noted that the factor k appears only in the explicit relation 
between x and wp», (57), and thus the derivation of equation (53) is not 
affected by the presence of k, although in this case w, in (29) must be given 
by (57). It is convenient to introduce a new variable x’ defined by 


7 a dae (59) 
— Rap2\(1—f)/48 
then (57) becomes 2’ = wi(3 ) : (60) 
Equation (60) gives an approximation to w in region (3) which ceases 
to be valid when x’ approaches 1. Putting 
2’—l=t (61) 


and assuming ¢ < 1, we get an approximate expression for the supersonic 
branch of (60) near t = 0: 


wy ~ 14+{—2(1—B)t}t (t < 0). (62) 
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Near x’ = 1, we can use an approximation to (24) already used by Wu (3) 
in the same region in the two-dimensional case. This approximation is also 
the same, in principle, as that used by Carrier (5) for finding an approximate 
solution to the relaxation oscillation problem. Using the transformation 


w = 1+atv*(n*), t = ain*, (63) 
equation (24) reduces to 
v*” + 2v*v*’+-2(1—8) = O(a). (64) 
Integrating and omitting the a‘ term, we find 
v®’ 4+ v¥2 4 2(1—B)n* = 0. (65) 
Putting v* = {2(1—f)}42, n* = {2(1—B)}-*n, (66) 
equation (65) becomes v’+v?+n = 0, (67) 


and its solution is of the form 
v = Z'(n)/Z(n), (68) 
where Z = (—£)*[K,(s)+ BK_,(—s)], 8 = 3(—7)!, (69) 
and A, is the modified Hankel function. Considering the boundary layer, 
where —€ > 1 and —t< 1, we can determine B by (62) which is valid 
in —t<1. For large values of —&, v~ —(—»)! if B+~O whereas 
v~ (—»n)' if B= 0 and 
Z = (—7)*K,(s) = (ant/v3)[ 4(s’) +-J_,(8’)] (70) 
with 3’ = $y}. 
This approximate solution, v*, in region (2) ceases to be valid when 7 
approaches ,, which is associated with the smallest zero of the function 
in the square brackets in (70). Near 7 = , 


a aah ee Sh. (71) 


, 


7—-N as | 
where 2} is 1+{2(1—B)}-*a*n,. 


v~ 





In the adjacent region (1), w varies rapidly from the value near 1 to 0 
in a narrow region and the approximate equation (45) 


etot” =C (72) 
can be used again, where é is defined by 


x’ = a, +ag. (73) 
The arbitrary constant c is determined by comparing (72) with (65) near 
x’ = 2x}; this gives 


c = 2—2(1—B)t, = 2—(ka)s, (74) 
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where ¢, = {2(1—)}-ta'y, and we put k* = {2(1—8)}'n,. With the value 
of c in (74), (72) is transformed into 


dw , w*—{2—(ka)}w+1 _ 





0, 
dé w 
which gives on integration 
E+c’ = (ka)-* tan-{(ka)-4(1—w)}—log(1—vw), (75) 


neglecting small quantities. Considering a region where (l1—w)-"at = 0(1), 
(75) reduces to 


E+c’ ~ (ka)-*(47)— —_— log(1—w). 


Comparing this with (71), which gives 1—w ~ —€-!, we get 


c’ = (ka)-*4n 
and & = (ka)-4{tan-{(ka)-*(1—w)}— }7]—log(1—vw). 
Using 2’ instead of £, we have 
a’ —a, = k-tai{tan-'{(ka)-*(1—w)}— }2]—a log(1—vw). (76) 
When w = 0, a’ ~ x, = 14+{2(1—B)}-ta'n. 


The typical solution terminates at a distance of order a from x’ = 1. This 
agrees with the result obtained by Wu (3) for the case of two-dimensional 
sink type flow. When a becomes zero the typical solution is as follows: 
w in 2, <x <2, (region 5), approaches w,(x) (subsonic branch) with 
jw—w,| = O(a); in |x»—ax| = O(a) (region 4), w approaches the straight 
line from wp; (= Wo(Xp)) tO Wg (Wo, Wog = 1). Inazy < x < kora,/k <2’ <1 
(region 3) w approaches w,(x’) (supersonic branch) with |jw—w,| = O(a); 
in |x’—1| = O(a) (regions (2), (1)) w approaches the straight line from 1 to 0. 

The typical solution for the case of a < 0 (sink type flow) can be found 
in the same way as above. Qualitatively, there is little difference between 
these two cases except for the fact that the roles of subsonic and supersonic 


branches in the inviscid solution are interchanged as we have seen in 
section 3. 


5. Other solutions 


When 2, in (40) approaches 1, the approximate solution (35) for f is no 
longer valid, since xo, x;, x2 in (33) become of smaller order than 1. Instead, 
we have an approximate equation (63) with (66), (68), (69), valid near 
x = 1. Since the subsonic branch of wy near x = 1 is expressed by 

Wy ~ 1—{—2(1—B)t}# where x—1 = t#, 
the arbitrary constant B in (69) must be chosen as non-zero in this case 





oil 
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and the solution curves near x = 1 become approximately as shown in 
Fig. 5, the exact shape depending on the value of B. 

For a more accurate determination of the values 
of B, we need a more precise approximation for w 
than w, near x = 1. In the case of two-dimensional 
sink type flow, Wu (3) used Lighthill’s method for !~; 
the purpose and this seems likely to be the case in 
the present problem also. 





For x > 1 the solutions approach a solution 
which is the continuation of the supersonic branch. 
At the limit of a +.0, this type of solution becomes 
the subsonic branch of the inviscid solution and the straight line. 

The special type of solution which corresponds to that in Fig. 3 (ii) is 
obtained if we take V2 A = —e-*@ instead of 


V2A = e-%', (40) 





1 > 
7,01 T 
Fic. 5 
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SUMMARY 


The steady flow of viscous incompressible fluid past an infinite flat quarter plate 
is considered, one edge of the plate being parallel to the main stream. Boundary 
layer equations are derived in a system of parabolic coordinates, and expansions in 
series provide a set of ordinary differential equations to determine the flow near 
this edge of the plate. An iterative procedure is described to solve these equations ; 
the approximations thus obtained indicate a decrease from the normal flat-plate 
value in the thickness of the boundary layer at the edge, and there is a singularity 
in the density of the skin friction at the edge. Both of these effects were predicted 
in previous investigations on a related problem. 


1. Introduction 


HirHERTOo the study of the effects of corners and edges on boundary layers 
has been concerned mainly with the Rayleigh problem. In this problem a 
viscous incompressible fluid bounded by two infinite half plates which 
intersect to form an edge (or corner) is set in motion by the sudden uniform 
motion of both plates parallel to their line of intersection; the determination 
of the velocity of the fluid depends on the appropriate solution of the heat 
conduction equation. Howarth (1) and Sowerby and Cooke (2, 3) have 
considered various cases, and have made use of a hypothesis due to Rayleigh 
to deduce qualitative information about the related steady flow problems. 
In this paper a direct attempt is made to consider certain features of the 
boundary layer for the steady flow of fluid past an infinite flat quarter plate, 
the main stream flow being parallel to one edge of the plate. The presence 
of an edge or corner leads to considerable complication in the form of the 
boundary layer equations, and the determination of the whole flow in the 
layer would be difficult. The results of the following investigation are there- 
fore restricted; we consider an approximation to the flow in the layer at 
the edge, in the plane of the plate, and hence determine the singularity in 
the skin friction at this edge. This singularity, and a certain decrease in the 
thickness of the layer at the edge, have been predicted by Howarth (1), and 
the results obtained here are in accordance with Howarth’s predictions. 
(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 3, 1958) 


eo | 
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2. Equations of motion in parabolic coordinates 

The flat quarter plate is taken to occupy the quarter planez = 0,2, y > 0, 
in a system of cartesian coordinates (x, y,z), the main stream velocity of the 
fluid being constant with cartesian components (U,0,0). It is convenient to 
replace y and z by parabolic coordinates a, 8, defined by the relations 


y=a®—f*, z= 208 (2.1) 


since the coordinate system (x,«,8) will later allow the formulation of 





z 
v , 
f Bb Bo 
. / 
te 
A=A>, 
a=0O / = 
B=o y 
Aa="-Xo 


Fic. 1. Diagram showing the system of coordinates ; the 

section is at right angles to the direction of the main 

stream. The coordinates &, 7 are given by the relations 

€ = a(U/vx)t and » = B(U/vx)t, so that in the diagram the 

curves a = const, 8 = const are also curves £ = const, 

7 = const, respectively 
boundary conditions in terms of a single coordinate. This system is an 
orthogonal system of coordinates, the surfaces a = const, 8 = const being 
parabolic cylinders as shown in Fig. 1. 
In this system of coordinates, let the components of the velocity v in the 

fluid be (u,v, w). The equation of motion for a viscous incompressible fluid 
in steady motion is 


} grad v?—vAcurlv = — : grad p—vcurl curl v, (2.2) 
p 


while the equation of continuity is 
div v = 0. 


In the coordinate system (x, «,8) this last equation becomes 


? (itu) 42 (hv) 4+ 4 (hw) = 0, 
Cx 


C 
Cx ep 
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where h, = 2(a?+?)!, (2.3) 


and hence, with the introduction of functions f (x, «, 8), ¥(x, «, 8), the general 
solution of the equation of continuity is 


htu=fp, hv=—ty bhw=da—-te (2.4) 
a suffix, as usual, denoting differentiation with respect to the variable 
indicated. The vector equation (2.2) then gives the following equations 
of motion: 


l ¢ n 
Py P(x tUpe+terg)» (2.5) 


p cx : ht 


] ] 
uu, — 7248 Uy + he Ug(ih —f,) — 
1 1 


1 l 4 ‘ 
— upg, “- i Pea— he (b.—Sf, bp = it (i+ (Ya—Se)*} 


Lé 8 | | 
— a {ga 38—Sra8) — Ta an +hpp—fua)+ Vice + Yea) 
(2.6) 
] 48, 
U(r ~Jan) = 52 Yaa Sra) 2 h? (ta —fr) (Yup —f2p)— hi 4 (at (ta —f,)"} 


1 
3 Op Sa ) 
“2 ep hs 2 (Yoana t $380 —S rac) — hi (oat gg Sead —Upet Pore Seae- 
2.7) 
3. The boundary layer equations and boundary conditions 
The above are the complete equations of motion, and we could proceed 
from these to deduce boundary layer equations in the usual way. How- 
ever, a solution of the problem in which y and z appear only in the combina- 
tion y/x' and z/x! may be expected; accordingly Carrier’s method (4) may 
be used to deduce the particular forms of the boundary layer equations 
which are required. That is, we now write 


s=(h eG 9 


(3.1) 


and 
f = (Uva )Hg(E, n) +-X9y(E, 0) + X*gelE, n)+..-} 
b = vi h(E, n)+X9d,(E, 7) +X*h0(E, n)+..-} (3.2) 
p = pU*{pl€, n) +X pil, n)+X*polE, n)+...} 


It is also convenient to write 


| 


u = Uiqé, n)+Xq(€, 9) + X*qe(E, n)+...}; (3.3) 





THE BOUNDARY LAYER FLOW ALONG AN EDGE 





so that, from (2.4), (3.1), and (3.2) we have 
4(£2+-n*)q = ~, ete. (3.4) 
U] 


If these expressions are substituted in the equations (2.5)-(2.7), and the 
coefficients of the various powers of X are then equated to zero, we obtain 
the following equations. From (2.6) and (2.7), 


Pe Pe  _. Mr _. 
ik hk tak = 0, (3.5) 
(€2+- ?)?q 7, (thet 1b) + (€?+ ”°)d, $én— 


—(E?+- 9?) ban(de—k) — Eth + (be—k)?} 
= (P+ 9 PP (E+ 1°) been + ban — ken) + 21 Geet+ Fun —Ke), (3-6) 
(é?4 vase (Ebe+ by) (e+ Eke tnky)| + (E+ 1b Beek) — 
— (E+ 9?)(pe—k)(hen—hy) + 0 b5+ (be—h)} 
= (7+ 9 PP —(E2+- 0?) (beeet- Senn — Kee) + 2E (Geet bnn—Ke), (3.7) 


where P = 4p,—q:—9, | (3.8) 
and 4k = 39g—£9:—n9, 
From (2.5), using the results of (3.5), we have 

— (E+ 9)E94¢+- In be— 39 + 1E9¢)—by Ve = eet Ione 3.9) 


There are further equations for g,, ¢,, etc., which will be neglected. The 
above equations, therefore, will describe the flow over the plate provided 
X is small; that is, at sufficiently large distances downstream from the 
leading edge of the plate. Equations (3.6)-(3.9), together with equation 
(3.4), are thus appropriate boundary layer equations for this problem. 
The expressions for the velocity components are, from relations (2.4), 


(3.1), and (3.4), 4(€2+ n*)u a Ug, (u =: Uq) 


vU 
wis ath tat ake (3.10) 


2+ thw = (=) eb 


It is clear that the above equations of motion are much more involved than 
the boundary layer equations for two-dimensional flow past a flat plate, 
even when allowance is made for the complications introduced by the choice 
of coordinate system. In particular, in the two-dimensional problem the 
component of velocity normal to the plate is obtained by simple integration 





rm ’ 
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of the equation of continuity; in this problem, however, the transverse 
components of velocity depend on the determination of ¢, which satisfies 
the fourth order partial differential equation obtained by eliminating P 
between equations (3.6) and (3.7). 

It remains to formulate boundary conditions for u, v, and w. The surface 
of the plate is defined comletely by the single condition » = 0; further, 
as 7 -> 00 (for any fixed . a... of 2) we move everywhere away from the 
plate into the main stream of the fluid. The surface = 0 is the plane 
extending away from the edge of the plate into the fluid, and coplanar 
with the plate. As € > -.0o we move away from this edge across the plate 
along its upper and lower sides respectively. This is illustrated in Fig. 1. 
The boundary conditions are therefore 





u=v=w=0, when 7 = 0, 
and u>U, asyn->o. 


Since ¢ satisfies a fourth order differential equation it should be possible 
also to specify boundary conditions on v and w as » > 00, namely 


v,w>0, asn>o. 


For the reason given above, it is not possible to impose this type of boundary 
condition in the two-dimensional problem. From the expressions (3.10), the 
boundary conditions will be satisfied provided 


q = ,, = de —k—0O when y= 0 | 


(3.11) 
q>1, $,,¢ 


z—k—>0O as n> \ 
$s 


4. The flow near the edge of the plate 

The feature of special interest in this problem is the flow in the vicinity 
of the edge of the plate, where € is small. We seek therefore a solution 
of the boundary layer equations in the form of a series of powers of &, the 
coefficients being functions of ». It is clear that g, g, and P must be even 
functions, and ¢ an odd function, of €. 


At this stage it is convenient to introduce a change of scale in the 
coordinate system by writing 


cf = g, | 
cp=m I 
where the constant ¢ will be determined at a later stage. Hence, if P is 
eliminated between the equations (3.6) and (3.7), and we write 
q = 3C*{Mo( 71) + £5 Mo( m1) +E} Ua) +-.--} 
_ eas §e{ fol) + £3 folm) +t falm) +---} ’ (4.2) 
¢= £1 41(1) +E) (71) +83 (1) +--- 


(4.1) 
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we obtain a set of equations to determine the functions which appear as 


coefficients of powers of £,. The conditions (3.11) give the following 
boundary conditions on these functions 





U=U—u=..=—0 ) 
= %,=—...=0 when 7, = 0 
ti—fo = Y3—bfe =... = 0 
/ (4.3) 
Uy > 3/c4, Uy, Ug... > 0 
Wh, Wg)... > 0 as 7, > 0 
t,—fotinto t3+i(mfo—fe) > 0 Ee 


and we may complete the definition of g by taking f,(0) = 0 (r = 0, 2,...). 
The first pair of equations is 


Uo+(fo—Yy)Mo+ 2u, = 0, (4.4) 

fo = 1%» (4.5) 

and although uw) may be eliminated easily to provide a single differential 

equation it is preferable, for numerical calculation, to retain the equations 

as they stand. The remaining equations become progressively more 
complicated. 

It is clear that the process of solving such equations is complicated by 
the appearance of higher order terms in the equations. For example, 
equation (4.4) contains the terms ys, u, and 2u,. The proposed method of 
solution, therefore, amounts to an iterative procedure, in which the expan- 
sions (4.2) are first developed in stages. The higher order terms are thus 
neglected in the first place in the equations, which are then solved to give 
approximate values of the various functions; these approximate values are 


then used in the complete equations to give second approximations for the 
functions. Thus we begin with 


q = 3¢*uo(m) 


g = s¢folm) |}; (4.6) 

$= 0 
so that the terms which are independent of £, in equations (3.4) and (3.9) are 
zero provided urturf, = 0, (4.7) 
fo = ni Mo (4.8) 


and the boundary conditions are 


ug(0) =0, fol) = 0) 


(4.9) 
Uy > 3/c* as n> 


The reason for the change of scale in the coordinates £,7 is now clear. 





—————— 
‘ 
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Following Tépfer’s method (5) for the solution of the Blasius equation, 


the numerical solution of the above equations may be simplified by replacing 
the last boundary condition by the condition 


u,(0) = 1, (4.10) 


and the scale factor c is then determined by writing 


act = 3, (4.11) 

where a is the asymptotic value assumed by u, for large values of 7. 
Equations (4.7) and (4.8) were solved numerically by the Adams- 
Bashforth method, and the results of this first approximation for wu) and 


f, are shown in Table 1. This table may be extended by use of the asymp- 
totic solution of the ev" ations, which is found to be 


Uy = a—A [ exp(—art+-ba) dx 


™m 





. ‘ (4.12) 
fo = jani—b+}A | (a? — 8) exp(—,,ax*+ ba) dx 
™ 
in which the values of the constants are 
a = 1-68465, b = 0-8922, A = 0:345. (4.13) 


The relation (4.11) now gives c = 1-15519. 
With these first approximations for u, and f,, the next stage is to extend 
the expansion (4.6) by introducing y,, so that 
q = jc*uo(m), 9 = §¢folm), = &,%4(), (4.14) 


where wp, fy are now regarded as already determined. The coefficient of , 
in the equation for ¢ is then zero provided 


4\ mw, [s ; ” oop 
+ for —S)NE +m 2 fod) Hwa = 0, (4.18) 
1 


and the boundary conditions for ¥,, from (4.3) and (4.12), are now clearly 


$40) = ¥(0) = 0 


, . 4.16 
yi, > 0, yy, > —b = —0-8922 as n, > oo | ( , 


Equation (4.15) is, of course, non-linear in y,. But it should be remembered 
that this method can yield approximations only for the various functions. 
In fact the first approximations for u, and f, have been derived by neg- 
lecting #, in the coefficient Y,—f,; further, the full equation for y, (which 
involves the higher order terms) introduces further terms on the right- 
hand side of the equation. It is consistent with the iterative procedure, 
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TABLE 1 
Values of the functions Ug, fy, 4, Ug, and f, 
(see equations 4.7, 4.8, 4.17, 4.19, 4.20, 4.4, 4.5) 
































Second First 
First approximation approximation approximation 
ui Ng M6 fo Ug So —H, —H 
° ° 1*00000 ° ° ° ° ° 

0°05 0°05000 | 100000 ° 00500 ° 0°0003 | O*O107 
o'10 O*10000 | 1*00000 | 0°00003 | 01000 ° oroorr | 0°0213 
Ors O*15000 | 100000 | 000013 | 0°1499 | O*000I | 0°0024 | 070320 
0°20 ©*20000 | 0°99998 | O:00040 | 0°1997 | 00004 | 0°0043 | 0°0427 
O25 0°25000 | 099995 | 000098 | 0°2494 | O°0009 | 0°0067 | 0°0534 
0°30 0°29999 | 0°99988 | 0°00203 | 0°2989 | O'00I19 | 00096 | 00641 
0°35 ©°34999 | 0°99974 | 0°00375 | 0°3482 | 0°0036 | o-o13r | 0°0747 
O40 | 039997 | 9°99949 | 0°00640 | 0°3973 | O'0062 | o-or7I | 00855 
O45 0°44993 | 099908 | O°01025 | 0°4461 | O°0100 | 0-0216 | 00963 
o"50 0°49987 | 0°99844 | 0°01562 | 0°4945 | O'or52 | 00267 | OvI07I 
oss 0°54977 | 0°99749 | 0°02287 | 75426 | 0°0223 | 0°0323 | o-1180 
060 059961 | 099612 | 0°03239 | 0°5903 | 0°0316 | 0°0385 | 01289 
0°65 | 0°64937 | 0°99422 | 0°04461 | 0°6374 | 0°0436 | 0-0452 | O*1400 
0°70 069902 | 0°99163 | 0°05999 | 0°6840 | 0°0586 | 0-0525 | O1512 
O75 0°74852 | o-9g8821 | 0°07903 | 0°7300 | 00771 | 00604 | 0°1625 
080 0°79783 | 0°98376 | 010228 | 07752 | 00996 | 00688 | 0°1739 
0°85 0°84688 | 0°97808 | 013029 | 0°8196 | 0°1267 | 0°0777 | 01855 
o"g0 089561 | 0°97094 | 0°16367 | 08631 | 01588 | 00873 | 01974 
O95 0°94394 | 096210 | 0°20305 | 0-9g056 | 0°1966 | 00975 | 0°2095 
1°00 0°99179 | O°95131 | 0°24908 | O°9470 | 0°2406 | 0°1083 | 0°2218 
1°05 1°03904 | 0°93832 | 0°30246 | 0°9g872 | 02913 | O-11Q97 | 0°2342 
I'Io 108558 | 092285 | 0°36387 | 1°0260 | 0°3494 | 01317 | 0°2468 
Iris 1°13128 | 0°90466 | 0°43405 | 1°0634 | 0°4155 | 0°1443 | 0°2594 
1°20 1°17599 | O8835r | 051373 | 1°0993 | 04900 | 01576 | 02719 
1°25 1°21958 | o85g921 | 0°60365 | 1°1335 | 0°5738 | O1715 | 02843 
1°30 1°26186 | 0°83160 | 0°70454 | 1°1659 | 0°6671 | o-1861 | 02962 
1°35 1*30268 | o-80060 | 081714 | 1°1964 | 07708 | o-2012 | 0°3075 
1°40 1°34186 | 076619 | 0°94219 | 1°2250 | o-8852 | 02168 | 0-3181 
1°45 1°37924 | 0°72846 | 1°08038 | 1°2516 | 1-o10g | 0°2330 | 0°3276 
1°50 1°41466 | 0°68757 | 1°23240 | 1°2761 1°1483 | 02496 | 0°3359 
Iss 1°44795 | 0°64384 | 1°39889 | 1°2985 | 12980 | 0°2666 | 0°3429 
1°60 1°47900 | 0759766 | 1°58046 | 1°3187 | 1°4603 | 0°2839 | 0°3481 
1°65 1°50769 | 0°54957 | 1°77769 | 173369 | 1°6356 | O-3014 | 0°3515 
1°70 1°53393 | o'50019 | I°ggITO | 1°3530 | 1°8243 | O°3190 | 0°3529 
1°75 1°55769 | 0°45023 | 2°22115 | 1°3672 | 2°0266 | 0°3367 | 0°3522 
1°80 1°57896 | 0-40045 | 2°46827 | 1°3795 | 2°2430 | 0°3542 | 0°3493 
1°85 | 1°59775.| 0°35165 | 2°73285 | 1°3900 | 2°4736 | 0°3716 | 0°3441 
1*g0 1°61415 | 0°30460 | 3°01521 | 1°3988 | 2°7187 | 0°3886 | 03370 
1°95 1°62826 | 0°26003 | 3°31565 | 1°4062 | 2°9785 | 0°4053 | 0°3279 
2°00 164021 | 021858 | 3°63444 | 1°4123 | 3°2534 | 0°4214 | O°3171 
2°05 | 165017 | 018074 | 3°97181 | 1°4172 | 3°5434 | 0°4370 | 03050 
2°10 1°65835 | 0°14688 | 4°32800 | 1°4211 | 3°8489 | 0°4519 | O°2917 

2°15 1°66493 | O-11721 | 4°70322 | 1°4241 | 4°1701 | 0°4662 | 0°277 
2°20 1°67014 | 0°09174 | 5°09769 | 1°4264 | 4°5072 | 0°4797 | 0°2631 
2°25 1°67418 | 0°07038 | 5°51164 | 1°4282 | 4°8605 | 04925 | 0°2435 
2°30 1°67724 | 0°05285 | 5°94532 | 1°4295 | 5°2303 | 0°5045 | 0°2339 
2°35 1°67952 | 0°03882 | 639900 | 1°4304 | 5°6167 | 05159 | 0°2196 
2°40 1°68117 | 0°02786 | 6°87293 | 1°4311 | 6:0202 | 05265 | 02061 
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TABLE ! (continued) 















































Second First 
First approximation approximation approximation 
mh Up U5 fo Mo to hy — 
2°45 | 1°68235 | 001952 | 7°36745 | 1°4316 | 6°4410 | 0°5305 | O°1934 
2°50 1°68316 | 0:01333 | 7°88287 | 1°4319 | 6°8795 | 0°5458 | o-1815 
2°55 1°68371 | 000887 | 8°41955 1*4321 7°3360 | 0°5546 | o-1705 
2'60 1°68407 | 0°00574 | 897783 | 1°4322 | 7°8107 | 0°5629 | 00-1604 
2°65 1°68431 | 0700361 | g°55811 | 1°4323 | 83041 | O'5707 | O1513 
2°70 1°68445 | 0-00221 |10°16077 | 1°4323 | 88165 | 05780 | o-1429 
2°75 1°68453 | o-00131 | 10°78621 1°4324 | 9°3483 | 05850 | o-1351 
2°80 1°68458 | 0-00075 |11°43483 | 1°4324 | 9°8997 | O°5915 | o-128I 
2°85 1°4324 |10°4713 | 0°5978 | o-1218 
2°90 1°4324 |11°0632 | 0°6037 | O-1I59 
2°95 1°4324 [11-6759 | 0°6094 | OII05 
3°00 1°4324 |12°3098 | 0°6148 | 0°1057 
3°10 0°6249 | 00970 
3°20 0°6342 | 00896 
3°30 06428 00829 
3°40 06508 | 00768 
3°50 0°6582 | 00710 
3°69 06650 | 00658 
3°70 06714 | 00618 
3°80 0°6774 | 00584 
3°90 0°6831 00552 
4°00 0°6884 | 0°0524 
TABLE 1 (continued) 
First approximation 
ul Mg Sr uh ug fs ™ Mg fe 
° ° ° 1°0O | O°1008 | 0°5255 | 2°00 | 0°0317 | 20822 
0°05 | 0'0058 | o-0013 | 1°05 | 01032 | 0°5817 | 2°05 | 00262 | 21704 
O10 | O°OII5 | O'0050 | II | O'1050 | 0°6408 | 2°10 | o-0213 | 2°2582 
O15 | O'O172 | OrO11Z | 1°15 | 01062 | 07029 15 | or0169 | 2°3456 
0°20 | 00230 | 00200 | 1°20 | 0°1067 | 0°7680 | 2:20 | 00132 | 2°4325 
O25 | 00257 | 0°0314 | 1°25 | O-1065 | 08359 | 2°25 | ororor | 2°5190 
0°30 | 0°0344 | 0°0452 | 1°30 | O'1055 | 09066 | 2°30 | 0°0075 | 26050 
0°35 | O'O40I | 0°0617 | 1°35 | O°1037 | 0°9799 | 2°35 | 00055 | 2°6907 
0°40 | 0°0457 | 00807 | 1°40 | O*1012 | 1°0557 | 2°40 | 00039 | 2°77 
0°45 | O'O512 | O-1024 | 1°45 | 0°0978 | 1°1339 | 2°45 | 0°0027 | 28610 
O50 | 00567 | 0°1267 | 1°50 | 0°0937 | 12142 | 2°50 | o-0018 | 2-9459 
0°55 | or0621 | 01538 | 1°55 | 0-088 | 1°2965 | 2°55 | o-o0r2 | 3°0305 
0°60 | 0:0674 | 01836 | 1°60 | 0°0835 1°3804 | 2°60 | o:0008 | 3°1150 
0°65 | 00725 | o-2162 | 1°65 | 0°0776 | 1°4657 | 2°65 | o-0005 | 3°1994 
O°7O | 00774 | 0°2516 | I°7o | oO-0712 | 1°5523 | 2°70 | o-0003 | 372838 
O75 | or0B21 | o-2899 | 1°75 | 0°0645 | 1°6397 | 2°75 | o-0002 | 373681 
0-80 | 00866 | 03311 | 1°80 | 0°0577 | 1°7278 | 2°80 | o-0001 3°4524 
0°85 | o-0g07 | 0°3752 | 1°85 | o-o509 | 1°8163 | 2°85 | o-o001 
O90 | 00945 | 04223 | 1°90 | 0°0442 | I°g050 | 2-90 | 00000 
0°95 | 070978 | 04724 | 1°95 | 0°0378 | 1°9937 
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therefore, to neglect the non-linear terms in (4.15) in deriving a first 
approximation for %,. The equation then becomes 


_ 4 ” - 2 "19 rue 
py +(f— =i ¥ (i ni— foi int Ut, = 0. (4.17) 
1 1 


This equation, with boundary conditions (4.16), was solved by use of the 
finite difference form of the equation, with an interval of 0-05 in », (this 
method was used for the solution of all the remaining equations). The 
values of ¥, and y, are shown in Table 1; it will be seen that %, does not 
tend rapidly to its limiting value, and this may be expected since an 
asymptotic solution of equation (4.17) yields a complementary function 
which behaves like »~! for large values of »,. 

The next stage is the introduction of the functions u, and f,, for which 
we take 

q = je*(Uo+Fi Uy), 9 = §e(fo+ fife), $= fh, (4.18) 
cea Wat (fo dadua+ ni tot Wi lua + Juofe = 0, (4.19) 
fo = Ut 73 Uy, (4.20) 
with boundary conditions u,(0) = f,(0) = 0, u, > 0 as n, > 0. The caleu- 
lation showed that these boundary conditions are not sufficient to define 
a unique solution for u, and f,—the condition at infinity appeared to be 
automatically satisfied by any solution which satisfied the conditions at 
7, = 0. A unique solution was obtained by making uw, tend to zero expo- 
nentially, in the same way that wu, tends to the constant value a, and it is 
reasonable to suppose that wu, must also behave in this way. Table 1 shows 
the values of wu, and f,. 

The next stages in the calculation would be the determination of approxi- 
mations for W,, U4, fy, ete., but the functions already calculated are sufficient 
to provide a second approximation for u, and f,, which is the main purpose 
of this paper. We return now to equations (4.4) and (4.5), and solve them 
with the same boundary conditions as before, using for ys, and u, the approxi- 
mations which have been derived; the results are shown in Table 1. 

In this second approximation, the new values of the constants a, b, c are 


a = 1-4324, = 0-5819, c = 1-2030, (4.21) 


so that a shows a 15% decrease from its original value, 6 a 35°, decrease, 
and c a 4% increase. The change in 6 is substantial, so that the tabulated 
values of ys, are probably too large numerically. However, an initial change 
of this order in ys, (and also, probably, a similar decrease in u,) might well 
be expected, since we began by taking 4, = u, = 0 in the calculation of 
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the original value of a, where the change is now only 15%. The evidence 
suggests that further stages in the iteration would yield a slightly lower 
change in a than this. Table 1 shows that the percentage changes in u, 
and f,, for all values of y,, vary from zero to a maximum of 15. 

It is clear that it would be desirable to continue the calculations, both 
to improve the approximations and to obtain some information about the 
convergence of the iterative process, but the work becomes progressively 
heavier. For example, in order to complete the next cycle in the process, 
which is to obtain a third approximation for wu, and f, (and thus a new 
limiting value for y,), it would be necessary to solve six further differential 
equations of which three would be of the fourth order. This and further 
work, indeed, might well require the use of an electronic computer. As far 
as can be judged there is no difficulty in extending the calculations, other 
than the considerable labour which is involved. 


5. Discussion of results 

The above approximations give an indication of certain features of the 
boundary layer. If the coordinates 

y = (~\y, a= ("Re (5.1) 
vat vat 

are now introduced, then in the plane z = 0,y < 0(sothat £, = 0) u attains 
a value within 1°%, of its main stream value when , = 2°15 (first approxi- 
mation) and 7, = 2-1 (second approximation); the corresponding values of 
Y are —3-46 and —3-06 respectively. If we regard the 1°, criterion as 
defining the thickness of the boundary layer, then the corresponding value 
in the flow over an infinite flat half plate (the Blasius problem) is Z = 5-0, 
the coordinate z being measured along a normal to the plate. The ratio of 
the two thicknesses is therefore 0-69 for the first approximation, and 0-61 
for the second approximation. It is probable that the true value of the 
ratio lies between these values, and is close to 0-61. The results agree with 
Howarth’s prediction (1) that a certain decrease in the boundary layer 
thickness is to be expected in the vicinity of the edge. 

The behaviour of the skin friction on the plate near the edge is of chief 
interest here. If + denotes the component, measured in the direction of the 
main stream, of skin friction per unit area, then 





ou [U\aed., - 
~~ (=) will (=) pi Hol) +ePFug(0)-+...} 
““ [2=0 


—~,0[U\en Vu! (0) + 
sos tei (c) Gll+erug(o)+..} (KW >0), (5.2) 


where the prime denotes differentiation with respect to 7. This expression 
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indicates a singularity in the density of skin friction at the edge of the 
plate. The order of the singularity is that predicted by the use of Rayleigh’s 
hypothesis (Howarth (1)) which gives, as the leading terms in 7, 
_ W2 _/U\ATG), REL . 
r=it o(;,) (1+ 0-33773¥ +...) (5.3) 
where ['(3) = 0-9064. 


The limit of the ratio of the value given by result (5.2) to the value given 
by (5.3), as Y > +-0, is 0-7447 with the first approximation for c, and 0-9120 
with the second approximation. This last ratio is probably slightly higher 
than the true value. The second approximation also provides an estimate 
of the coefficient of Y* in (5.2); the value of u(0) is 0-1150, so that 

c2u,(0) = 0-1664. 


This is about one-half the corresponding coefficient in (5.3). 


REFERENCES 
. L. Howartn, Proc. Camb. Phil. Soc. 46 (1950) 127. 
. L. Sowersy, Phil. Mag. (7) 42 (1951) 176. 
——and J. C. Cooker, Quart. J. Mech. App. Math. 6 (1953) 50. 
G. F. Carrier, Quart. Appl. Math. 4 (1946) 367. 
C. Toprer, Z. f. Math. u. Phys. 60 (1912) 397. 


a PON 





en 


‘pac hse tari 


A SOLUTION OF THE NAVIER-STOKES EQUATIONS . Aa 

ILLUSTRATING THE RESPONSE OF A LAMINAR 

BOUNDARY LAYER TO A GIVEN CHANGE IN THE 
EXTERNAL STREAM VELOCITY 


By J. WATSON (National Physical Laboratory, Teddington, Middlesex) 
[Received 5 September 1957] 


SUMMARY 


Exact solutions of the Navier-Stokes equations are derived by a Laplace-transform 
technique for two-dimensional, incompressible flow past an infinite plane porous wall. 
It is assumed that the flow is independent of the distance parallel to the wall and 
that the velocity component normal to the wall is constant. A general formula is 
derived for the velocity distribution as a function of the given free-stream velocity 
and, from this, general formulae for the skin friction and displacement thickness are 
obtained. 

Among the particular cases considered are (i) flow in which the free-stream velocity 
is changed impulsively from one value to another, (ii) uniform acceleration of the 
main stream from a given value of the velocity, and (iii) flow in which the free-stream 
velocity consists of a mean value with a superimposed decaying oscillation. In the 
first case a secondary boundary layer of a Rayleigh impulsive flow type is created 
immediately after the impulsive increase in free-stream velocity and this secondary 
layer increases in thickness with time until the final steady state is reached. In the 
second case also, a secondary boundary layer is formed at the instant when the 
free-stream velocity begins to accelerate and ultimately the skin friction consists 
of a quasi-steady part and a part proportional to the product of the free-stream 
acceleration and a ‘virtual mass’ equal to the mass of fluid in the displacement area. 
In the final case the fluctuating skin friction is found to have a phase advance over 
the fluctuating free-stream velocity. For sufficiently large frequencies this phase 
advance is }7 and the amplitude of the skin friction fluctuation is proportional to 
the square root of the frequency and to the amplitude of the free-stream velocity 
fluctuation. 

























1. Introduction 


LiGHTHILL (1) has investigated the effect of fluctuations of the external 
stream on the boundary layer of a two-dimensional body. In addition, 
Stuart (2) has obtained an exact solution of the Navier-Stokes equations 
for the case of fluctuating flow past an infinite plane wall with constant 
suction, and has compared some of his results with Lighthill’s. This paper 
is an extension of Stuart’s work to the case of generai unsteady flows past 
the wall, but with the suction still held constant. 

The generalization of the velocity field is effected by means of a Laplace- 
transform technique; many forms of the time-dependent external stream 
[Quart. Journ, Mech. and Applied Math., Vol. XI, Pt. 3, 1958] 
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velocity can be considered, since the mathematical conditions which the 
velocity must satisfy are only very weak. General formulae are given for 
the boundary-layer velocity distribution and for the skin friction. Among 
the particular cases considered are (i) flows in which the external stream 
velocity goes from one value to another, either impulsively or with accelera- 
tion, (ii) the case of a uniform acceleration from a given value of the velocity, 
and (iii) the case of an external stream which consists of a mean value with 
a superimposed decaying oscillation. 

Immediately after an impulsive increase of velocity at time t = 0, a 
secondary boundary layer is created whose thickness is of order ,/(vt), where 
v is the kinematic viscosity. Outside this thin layer the main boundary 
layer behaves as if it were inviscid, and at each point has its initial value 
plus the same impulse as is given to the external flow. As time progresses 
the secondary layer increases in thickness and gradually distorts the flow 
until the ultimate steady-state boundary layer is reached. It follows that 
because of the presence of the very thin secondary layer, the skin friction 
at first rises to a large value, and then decays to its ultimate steady value 
as the secondary layer increases in thickness. The initial thin secondary 
layer is of the form of a Rayleigh impulsive flow, and it can be inferred that 
when any boundary layer is given an impulsive shift from one velocity to 
another, the flow will react as though it were inviscid except within a thin 
secondary layer. In a separating boundary-layer flow it can further be 
inferred that the large rise in skin friction would have an appreciable effect 
on the transient position of the point of separation. 

In the case of a boundary layer with an external stream accelerating 
uniformly from a given value at t = 0, a secondary boundary layer is formed 
at small times whose thickness is of order ,/(vt), and for a free stream 
U,(1+-kt) the skin friction initially grows like t!. This result follows because 
outside the thin secondary layer the fluid behaves as though it were inviscid 
and has its velocity increased by U, kt; thus the increase in velocity rises 
from zero at the wall to U, kt in a distance of order , (vt), leading to the 
above skin-friction growth. The skin friction at large times is found to 
consist of two parts, the first of which is the same as if the skin friction 
were steady at the instantaneous free-stream speed. Together with this 
quasi-steady skin friction is a part which is proportional to the product 
of the free-stream acceleration and a ‘virtual mass’ equal to the mass of 
fluid in the displacement area; this second part of the skin friction arises 
because, through the action of viscosity, the pressure gradient does not 
accelerate the fluid in the boundary layer as much as the free-stream and 
consequently there is a virtual acceleration in the negative direction, 
corresponding to a drag. By combining two accelerations at different 
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times, the case of a fluid accelerating from one velocity to another is 
considered. 

The final case dealt with is that in which the free-stream velocity is 
constant for ¢ < 0 and has a superimposed decaying oscillation for t > 0. 
The skin friction is found to have a phase advance over the velocity fluc- 
tuation. For sufficiently large frequencies or, in other words, for a small 
logarithmic decrement, the skin friction behaves as though the free-stream 
oscillation were a non-decaying one with its instantaneous amplitude. Then 
the amplitude of the skin friction fluctuation is proportional to the product 
of the square root of the frequency and the instantaneous amplitude of the 
free-stream velocity fluctuation; moreover, the phase advance is }7. In 
this quasi-harmonic state, conditions are analogous to those for pure oscil- 
latory flows of high frequency (Lighthill (1)). 

A final point is that the velocity distributions obtained by E. J. Watson (3) 
for the flow past an infinite flat plate without suction are particular cases 
of the general solution given in this paper. 


2. List of symbols 


x coordinate parallel to wall. 
y coordinate normal to wall. 
t time. 
t, characteristic time. 
p pressure. 
u velocity in x-direction. 
v_ velocity in y-direction. 
free-stream velocity. 
U, constant reference velocity. 
v, (constant) velocity component normal to wall. 
p density of fluid. 
viscosity of fluid. 
v kinematic viscosity of fluid (y/p). 





{o(y) defined by equation (3.4). 
f(t), gly, t) defined by equation (3.5). 


T’ non-dimensional time (tv®,/4v). 
7, non-dimensional characteristic time. 
7 non-dimensional distance normal to wall (y|v,,|/v). 
© non-dimensional frequency. 
519 steady displacement thickness (v/|v,,|). ' 


6, displacement thickness, defined by equation (3.13). 


Tw 


skin friction. 
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rE 
erfx error function (; | e-=* ix) 
7 
0 
erfex l—erfz. 
_ zx d . 
C(x) Fresnel integral - aed, 
(27)* xt 
= 0 
zx 





S(x) Fresnel integral = ain z de ‘ 
(27)! at 
rs 0 








3. General theory 


Let x denote the distance along a two-dimensional, infinite, plane porous 
wall, y the distance normal to it, uw and v the corresponding components 
of velocity, t the time, p the pressure, p the density, and v the kinematic 
viscosity. For two-dimensional incompressible flow which is independent 
of x, the Navier-Stokes equations and the equation of continuity are 





a 2. 
a. ae 
et oy p\ex oy? 
+ = — (2) ; (3.1) 
ot p\cy 
——ia® 
oy } 


subject to the boundary conditions u = 0 at y = 0 and u > U(t)asy >, 
U(t) being the free-stream velocity. In order to obtain a steady solution of 
the boundary-layer type it is known that v must be a negative, non-zero 
constant (v,,). In the unsteady case also we shall make this restriction 
(Stuart (2)); then p is independent of y, —(1/p) @p/éx equals dU/dt, and 
the first equation (3.1) becomes 


iat iat (32) 
Stuart (2) considered a periodic velocity field of the form 
U(t) = U[1+ee] 
u = UlSo(y)+eli ye] ?, (3.3) 
v= 0, 


w 


where « is a constant, and found that 
f(y) = 1—exp(—y|v,,|/v) 
and f(y) = 1—exp(—Ay|v,,|/v), (3.4) 
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where h = }{1+.,/(1+4iwv/v?,)}. He also pointed out the existence of an 
exact solution of (3.1) of the formt 
) = Uf1+f(0] 
dia Ul Soy) +-9(y, t)] ? (3.5) 
v= v, 
where {,(y) is given by (3.4) and where f(t) is a given arbitrary function 
and g(y,t) is a function to be determined. The present paper deals with 
an investigation of this generalized case. 
Let us first make the transformation 
T = tv?,/4v = vt/482,, n= y|v,|/¥, (3.6) 
where 8,, is the steady displacement thickness (v/|v,,|). Substituting (3.5) in 
(3.2) and utilizing the fact that f, is given by (3.4) we obtain 
ég(n,T) _ 409: ,T) = f(T) + 42 *9(n, T) (3.7) 
oT on 7” 
with the boundary conditions g = 0 at y = 0 andg > f(T’) as yn > ©. 
In solving (3.7) we use the two-sided Laplace transform, which eliminates 


the independent variable 7’. It is defined by van der Pol and Bremmer (6) 
(with a change of notation) as 


i. 





#(p) = p | ePTy(T) dT. 
The inverse is AT) = iF | evT=(P) dp. 
a7 Pp 


A formal solution of (3.7) by the Laplace transform technique will now be 
derived. To do this we note that the transform of the first term on the 
left-hand side of (3.7) is 


x « 


p J: ror.) dT = =p fe e-?Tg dT = pG(n,p); 


. 
x —@ 


provided that e~?79(n, 7’) > 0 as T' > +00, where (7, p) is the gh 
of g. Similarly, the transform of f’(7') is pf(p), provided that e-?7f(7') > 0 
as T'-> +00. The Laplace transform of (3.7) is therefore 


dg dg oe eS 


8 
dy? dy 49 a4) 


+ After this work was completed, the author found that the solution when U(t) = 0 
for t < 0 has been given by Hasimoto (4). 

¢ There is also a solution of (3.1) in which the steady part of the free stream has a 
uniform shear. The unsteady part g of the velocity profile u is independent of the shear 
(Curle (5)) and so is given by the solution of this paper for zero shear. 
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with the boundary conditions 


Sib Conia. (3.9) 
g +f as 7 > 0 j 

If k = 4{1+-,/(1+-p)} and that branch of the square root is chosen which 
makes ,/(1+-p) = +1 when p = 0, then the function {(1—e-*") satisfies 
(3.8) and the boundary conditions (3.9); the real part of k is positive for 
all values of p with this choice of the square root. It follows that 


g = f(1—e-*”) (3.10) 
is the unique solution of (3.8) subject to (3.9). 
It is now necessary to determine g. We first note that (3.10) is the 
product of two functions of p one of which, f(p), has a known inverse, f(7’). % 
Thus the function g is best obtained by using the composition-product rule, 
namely that the inverse of 


@ 


t. " . 
pii(P)t(P) is | 2, (A)a,(7'—A) da. 


We identify f with z,, f with x,, and p(1—e-*") with Z,. We now need the 
inverse of the transform p(l—e-*"). Using two standard results we find 
that this inverse is given by 


c+ia 


4(7')— e~ jn) [ ePTe-iw+P) dp 
2m . 
c—ia . 
e*+ia F 
at J Pp 
Fa d * 2 
aS 3(T)—+ e-ine-T BW | ePAT in ne—Vv ©P (: =— c ) 
7? 2m p 4 
c’—ia 
/16T 3 
= 8(T)—S,e-e- Ha) nd. 
7? 7 iT! 8 


eth a ) 
or) — el 


e—Nne-Te- m16T 
where 6(7') is the delta function and H(T) is the unit function. The 
composition-product rule then gives 


nH (4A/n?) 
‘4ndi 


g(n, T) = fa (T— [aa —™ e-IneNe-n nail dn 





=4-"? [ -yer—aexp|—(r4+%)} a. (3.11) 
0 
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It is assumed that A (7) is bounded in and integrable over any finite 
range and the integral | A-if(T'—A)e~ dA (A > 0) is absolutely convergent 


for finite 7’. In such a case, it can be shown that g(», 7’) as given by (3.11) 
satisfies (3.7) and, at points where f(7') is continuous, that it satisfies the 
boundary conditions. Consequently (3.5), where g(y, t) is given by (3.11) and 
(3.6), satisfies (3.1) and the corresponding boundary conditions at points 
of continuity of f(7'); on physical grounds the solution is unique. Functions, 
f(T), which do not satisfy the above conditions may, however, lead to a 
solution; in such cases it must be checked that (3.11) satisfies (3.7) and the 
boundary conditions. From equations (3.4) and (3.11) the velocity distri- 
bution, (3.5), is given by 
U(T)|U, = 14+f(T) 


ne“ 


at [ a-yer—ayexo{ — +74) da }, 


u/U, = 1—e+f(T)— 


0 
y= ¥ 
(3.12) 
where 7’, » are related to t, y by (3.6). It can be shown that if the free- 
stream velocity ultimately becomes steady, then the velocity distribution 
reaches a form similar to the undisturbed velocity distribution; that is if 
lim U(T) = U(#), where U(0) is a finite constant, then 


To 
lim u(y, 7’) = fo(n)U (2). 
To 


Similarly if lim U(T) = U( oo then 
T—>-—« 
lim u( T) = Lo(n)U(—c0) 
under similar conditions. These results are proved in Appendix A. Par- 


ticular cases of the general velocity distribution are discussed in the 
following sections. 


The displacement thickness, 5,, may be defined by 


x 


J 1 | dn, (3.13) 


which will now be evaluated ink’ (3.12) under the assumption that the 
integral 


@ 


{ A-4f(T—A)e—> dA 
i) 
is absolutely convergent for finite 7’. From (3.13) 


ee u(n, =| a " 





U, U, 
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which yields, after a little algebra and use of (3.12), 


1 9 i 9) @ 
eo Beene -t —))e- dr —5 oan 4 
8, |v|/v cage ta | f(T—A)e-> ar 2 fr Ayerfed a 
0 
(3.14) 
It now remains to derive a formula for the skin friction, 7,,, which is 
given by 
oiiliie) ws je.) /en| = lim “9 7) (3.15) 
Uy n=0 770 Uy 


since u—> 0 as 7 > 90. From (3.12) and by making use of equation (B,2), 
this becomes 


rl PUalMel = VLD) +75 | MM(T)-S(T—Aje* A, (3.16) 
0 


provided that the integral in (3.16) is absolutely convergent. This is so if 
at the point 7’ 


f(T)—f(T—aA)| < C\A\¥ for |A| <3, (3.17) 
where C is a positive constant and k > 3. We note that a function which 
satisfies (3.17) is necessarily continuous at the point 7’; while a function 


which has a finite derivative at 7’ satisfies (3.17) and the skin friction is 
given by (3.16). 


4. Applications of general theory 


In the following examples all the conditions on f(7') = U(T’)/U,—1 are 
satisfied in order that the velocity distribution, the displacement thickness 
and the skin friction be given by (3.12), (3.14), and (3.16) respectively. 


4.1. Periodic velocity field 
Application is made to the case f(7') = ee7, where Q = 4yw/v?,, in 


order to demonstrate that certain results given by Stuart (2) can be derived 
from the results given in section 3. The second equation in (3.12) becomes 


3) 


u : np G78? heal { 7? \ 
— ]—e-1-+ eink _-iaT ie-iMexp! — A+—.}} da 
U, ae ss 473 [> _—s ( 5) 


wo 
-4 
—= ]—e-14 ei Tegan 1° ° N-He—(a*A+v2/~) dy 
4nt : 
0 


where a = (1+iQ)! with re(a) > 0 and b = }». Using the evaluation of 
the integral given in Appendix B, equation (B, 4), we obtain 
yn 4rrt net , 
in - i To! p41 +inyt 
u/U, == Qing 74 e¢iMT — ¢gif ite e~ 8m +iQ) 


= 1—e-7+ eeiT(1 —e-"), (4.1) 
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where h = 4{1+-(1+iQ)*} and Q = 4vw/v?,, which is a result given by 
Stuart. 
From (3.16) the skin friction is given by 


Tp|pUs|%p| = 1+ee@? - +5 =¥ A-#(1—e-#e— Ad. 


0 


The integral in this equation may be integrated by parts to give 





A-He A_(1 Li Q)e (1+i0)A} dX 





T \ay - jl -ptQT mi 
Tw ply Uw =1 €€ a 


= Jt ¢eiNT —eei Tf] t+—4(1-+4 iQ) 
= ltee®?h, (4.2) 
where h = 3{1+-(1+iQ)}!}, by use of equation (B, 6) of Appendix B. This 


agrees with Stuart’s result. 


4.2. Impulsive velocity fields 

(a) Single impulse 

The case of a single impulse at 7’ = 0, namely 

U(T) = U, for T < 0; U(T) = U,1i+A) for T >‘ 
where A is a constant, is now considered. This corresponds to 
f(T) = AH(T), 

where A is constant and H(7’) is the unit function. The velocity is given 
by (3.12) and, because of the factor H(7'—A), the integrand is zero for 
A > T. Consequently it is only necessary to integrate from 0 to 7’ when 
T’ is positive; for negative values of 7' the integrand is identically zero. 
Thus, provided 7 + 0 the velocity has the form 


T 
u | Ane-1H(T) [a 
= ]|—e-7+- A(T) — AOE coat. Bs -4 we : 
7 é (T) it — | Mtexp| (+75) dn 
0 
273 ne in 


= 1—e~4+ AH(T)—AH(T) = elnx 
4h 


‘ u] 
x le nerfe( 7 Ti 7") 4 erfo( 7+ r")) 


from equation (B, 2) of Appendix B. Hence 


a= Jin e~Git 7 va : 7 
1—e-1 an(r)| je-nerfe( 7, 1) — jerfe( 75 +7")| (4.3) 


U 


0 


0 


for 7’ ¢ 0. Itis easily seen from this that, as 7’ > oo, u/U, > (l—e-7)(1+A) 
as we would expect from Appendix A with U(oo) = U,(1+-A). 
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The displacement thickness, 5,, as defined by (3.13), is given by (3.14) for 
T #0. Thus 


T 
, 1 ao e~ od 
. -o l 9 > 4 
5, /0,,|/0 ‘TT AMT)} [1+ is —2AH(7) { exfer a (4.4) 
0 
From Appendix B, equation (B, 5) gives 
T 
B awl 
| = = rierf T}. (4.5) 


« 


0 
Together with this result, a double integration by parts of the second 
integral leads to 
- 
[ erfe At dA = Terfe T!—(T/n4)e-? +} erf T! (4.6) 
0 
so that (4.4) becomes 
il |/y = a 
s {14+AH(T)} 
for T + 0. 
Finally, the skin friction, 7,,, given by (3.16) is, for 7’ 4 0 


1+AH(T){(27!/7)e-? +erf T!—2T erfe T*}| (4.7) 


| 
Te ply Cy| = 1+ AH(T) 24 [ e ‘dn 











4nt | Dl 
; 
,, , AH(T)| 2e-7 fs eA dy 
os 24 anc | 08 So. 
1+AH(T) + a nl? 7 : 
T 
oe y 
-X —A e o-A 
But “= a i* a : — = mt—nierf Tt (4.8) 
T 0 0 
from (4.5). Hence 
» 
e jell ben = 14 aH) +222 On n+ 2nterf 7%] 
4a 
T 
2 1p iter + = (4.9) 


for T 4 0. 
A graph of the skin friction, 7,,, for A = 0-2is given in Fig. 1. Eventually 
7,, reaches the steady value +, which, from (4.9), is given by 


t,[pUs|v,| = 1+-A. (4.10) 
It can be seen from Fig. 1 that at 7’ = 0 the skin friction initially goes to 


plus infinity and then decreases to its limiting value as T increases. The 
reason for these properties of the skin friction can be seen as follows. Just 
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before the main stream velocity is changed impulsively from U, to Uj(1+-A) 
the velocity distribution in the boundary layer is Uj(1—e-"). Just after 
the application of the impulse, that is for very small values of the time, 
it can be seen from (4.3) that the velocity distribution is of the form 


u/Uy = 1—e-7+ A—43A(1+-e-) erfe(n/47") 


which depends on the parameter 7/7. From (3.6) this parameter corre- 
sponds to'the parameter y/(vt)!. For values of y such that y > (vt)! the 
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Fic. 1. Behaviour of skin friction in impulsive flow 


velocity distribution is U,(1—e-”+-A), whereas for values of y of order 
(vt)', that is when x is very small, the velocity varies from zero velocity at 
the wall to the approximate value U,A in a distance y of order (vt). In 
physical terms, for very small times there is a secondary boundary layer 
on the wall, whose thickness is of order (vt)!. As t > 0 the same velocity 
change, U, A, must take place in an ever decreasing distance; thus the skin 
friction becomes infinite as t > 0 like 1/t!. As time progresses the secondary 
boundary layer increases in thickness and distorts the flow in such a way 
that the initial velocity U,(1—e-"+-A) approaches more and more closely 
the ultimate form U,(1+A)(1—e-"). A further fact worth noting is that 
the viscous effects of the impulse are confined to the secondary boundary 
layer; so that the flow outside the secondary layer may be regarded as 
inviscid as far as the impulse is concernedt and the velocity there increased, 


+ A similar phenomenon is noted by Lighthill (7) in his work on shox boundary- 
layer interaction. He points out that it is only in ‘an “inner viscous sublayer’ ... that the 
disturbances to the viscous forces are comparable with the disturbances to the inertial 
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for T’ > 0, by U,A, the impulsive increase in the free-stream velocity. 
Consequently, by letting 7'-> 0, we see that, for » > 0, the velocity is 
impulsively increased by U, A at 7’ = 0. It should be noted that, for small 
values of » and 7’, the secondary boundary layer is of the form 


u = U,Aerf (7/47), 


and is thus like the Rayleigh boundary layer for impulsive flow from rest. 

It is of some interest to define a non-dimensional time, 7., which is 
characteristic of the time taken for the skin friction to change from one 
steady value to another. In the case under consideration this can be defined 
as follows. Let 7, be the undisturbed skin friction, namely, 


7, = pU,|v,,|. (4.11) 
Then the characteristic time 7, is defined by 
T, = | (2=—Jar, (4.12) 
J \Te— 71 


where 7, is given by (4.10). From equation (4.9) this is 


r= | Visertms © \_ilar. 
. 2| (xT)! 
0 
On performing the integration we obtain 
T, = }. (4.13) 
From (3.6) this gives t, = 8%5/v, (4.14) 


and for a boundary layer of displacement thickness 0-1 cm this is t, => 0-07 
sec. The time taken for the flow to settle down to its ultimate state will 
be about three times this. It can be seen that the characteristic time, in a 
sense, is analogous to a displacement thickness. 


(b) Multiple impulse 

From section 4.2 (a) it is a simple matter to derive the velocity distri- 
bution for the case when the free-stream velocity is given a series of 
impulses. In particular consider the case 


U (T < 0) 
U(T) = (U(1+4A) (0<T<T) 
Uy (T) < T) 


forces’. Hence outside this layer, though still in the boundary layer, the disturbances to 
the viscous forces are negligible, and the interaction is quasi-inviscid. 
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so that 0 (1 <0) 
{(T)/={A (0<T<T% 
0 (T< T) 
* f(T) = f*(T)—f*(T—-TN), (4.15) 
where f*(T) = AH(T). 


From (3.12) it is seen that the part of the velocity distribution which depends 
on T is linear in f so that by (4.15) it is composed of the part arising from 
f *(7) together with the part arising from —f *(7’—7). But, from (4.3), the 
part arising from f *(7’) is 
$(n, T) = an(T)|1 2 Je-nerfe( ne m)— \ erfe( hy 4 r")| (4.16) 
for T #0. Similarly the part arising from —f*(7’'—7T) is —¢(n, T—T) 
for T ~ T,. Thus the velocity distribution is given by 
u/U, = 1—e-7+ 4(n, T)—¢(n, T—T) (4.17) 
for 7 ~0 and T +7, and where ¢(n, 7’) is defined by (4.16). Since 
d(, 7’) > A(1—e-”) as T' + oo then, from (4.17), u/U, > (1—e-”) as T' > 00 
| U(o) = Uy in Appendix A]. By similar reasoning the skin friction, + 


w? is 

given by Fpl PUp |p| = 1+0(7)—Y¥(T—T) (4.18) 
for 7 ~ 0 and T + 7%, and where (7) is defined by 
_ AH(T) e-? * 

wr) = ‘rere T+ (4.19) 


Finally, the displacement thickness, 5,, is given by 
8, Moll” = (1+ x(7)—y(T—T [1+ AH(7)—AH(T—T))] (4.20) 
for 7 A 0 and T + JT, and where y(7’) is defined by 


2T%e-T 
x(T) = an(n)|° ~~ + erf T!—27 erfe 1). (4.21) 
7 


A graph of the skin friction, 7,, for A = 0-2 and 7, = 0-04 is given in 
Fig. 2. For T < 1, ¥(7—T,) = 0 from (4.19), and (4.18) and (4.9) are 
identical so that, for 7’ < 0-04, Fig. 2 represents the same skin friction as 
Fig. 1. By the same reasoning as in section 4.2 (a), a decrease in free-stream 
velocity from U,(1+-A) to U, at 7 = T, causes the skin friction to go to 
minus infinity and then to increase to its limiting value. From (4.18) and 
(4.19) this steady limiting value, 7,, is 


7,|pU,|v,,| = 1. (4.22) 
Again, as in section 4.2 (a), a new secondary boundary layer is created at 
T = T, the thickness of which is of order v4(7'—7,)t for small positive 
(7’—T)) and as time progresses it increases in thickness and modifies the 
flow to the ultimate steady velocity U,(1—e-"), By the same reasoning as 
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in section 4.2 (a), the flow outside this new boundary layer may be regarded 
as inviscid as far as the impulse at 7’ = 7, is concerned and for 7’ > 7), the 
velocity there is decreased by U, A so that at 7’ = 7, the velocity for 7 > 0 
is impulsively decreased by this amount. 


2-0; 
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Fic. 2. Behaviour of skin friction in flow with two impulsive changes of velocity 


4.3. Accelerated velocity fields 

(a) Single acceleration 

Consider now an accelerated free stream with velocity of the form 

U(T)=U, forT < 0, U(T) = UJ1+KT] for T >0 
where &K is a constant. This corresponds to 
f(T) = KT H(T), (4.23) 

where (7') is the unit function. Note that, although f(7') is not differen- 
tiable at 7’ = 0, it satisfies (3.17) there. ae (3.12) and (4.23) yield 


ries 2 14-KTH(T)— a (T) fx (7 — dNexp! — (+77)] dn 


U, | 16A 


| 





T 
r , —} 
1—e-n- KT A(T) — A Ae , [ 2- texp| — Bi dd + 
3 


4 
KH(T)ne™ (4 f(y, 7° \\ 
0 


a 
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But from Appendix B, equation (B, 2) gives 


T 
2n inl e- u} 7 
| A- texp{ — (a+ a} dA = =melnle nerfe(t5 — 7") +-erfo( 7+ r")) 
0 


and equation (B, 1) gives 


T 

n\ e- . _ 
| A- texp{— (a+ is} dA = 4ret fe nerfo( 7, — 7) —exfo( 7, + )| 
0 


so that (4.24) may be written in the form 


— sa 7 7 \ 
= 1—e-7+ 4K H(7)| sr — ar rerfo( 7, — 7") + erfe( 15 + 7")} + 


tafe nerfe( 74-1") orf nt)|| (4.25) 


for all finite 7. It is easily seen that, for 7 > 0, as T’ > 0, u > U,(1—e-") 
and as 7’ +>, 


Uy 


u~ Ujfl—e-+ K[ T(1—e-) +} ne-"]}. 
From (3.14) and (4.23) the displacement thickness is 
; 1 2K A(T) 
\v..|/v = ——__—__. |] 4. * | -4(T’—A)e- dA— 
8; Vwp|/¥ tr + po fr ( ye 
0 


T 


—2K H(T) { (T—A)erfer a| 
0 


T 

1 2KT H(T) 

er F poems Bch eH —t eA oor 
aT ETAT +a is ings 


_ 2K =) 





r 

fm -\d\ —2KT H(T) | erfe At dA + 
0 
T 


49K r) | Aerfo At a|. (4.26) 


Of the integrals on the right-hand side of (4.26), the first is given by (4.5), the 
second can be derived from (4.5) and is 


T 
| Ate-A dA = Anterf T!—T'e-T, (4.27) 
0 
the third is given by (4.6), and the fourth can be obtained by a double 
integration by parts which, after using (4.27), gives 
T 
- 3 te-T 
[ AerfoAt dA = 47Perfe TH+ gerf TH 7 " _3Me" (4 98) 


4 27! 4nt 
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Substituting in (4.26) and simplifying, we obtain 








Y ae 1 4 
8, |Y%|/¥ = cere +e marx 
ip-T te-T 
(Te? Te morte TH4 Perf Ti—jerf TH] (4.29) 
| or 2r3 J 
for all 7’. 


Finally, from (3.16) the skin friction, 7,,, is, for finite 7’, 


rol PUe|te| = 14K THT) 45D) 


7 


T a) 
| | AeA ALT | A-le-A al (4.30) 
0 T 
But 


i ) 


[ Ate dd = [—20-e AP —2 [ Ate dA = 2e-P Tt 2at-+ Dat erf TH, 
T T 


by (4.8), while the first integral in (4.30) is given by (4.5). Thus 
Ty|pUy\Cy| = 14+ K A(T) 47+ ferf Tt+407%Tte-7+47 erf Tt] (4.31) 


for all finite 7’. A graph ofr,,/pU,/v,,| isgivenfor K = 4(Fig.3). Eventually 
7, reaches the form which, from (4.31), is given by 


74o/PUo|%o| = 1+K(T+}). (4.32) 
In Fig. 3 the gradient of the skin friction curve is seen to be infinite at 
T = Oand then it steadily decreases to a limiting value as the time increases. 
The reason for this behaviour of the skin friction at 7’ = 0 is similar to the 
reason for infinite skin friction in the impulse case (section 4.2 (a)). In this 
case, for small positive 7’, the velocity changes from zero at the wall to 
U{(l1—e-") + KT} in a distance » of order T*; so that 7,,/pU/v,,| = 14+CT* 
where C is a constant. That is, the skin friction grows parabolically for 
small positive 7’. As time progresses the secondary boundary layer, with 
thickness of order 7, continues to grow and, as it does so, the skin friction 
changes. Ultimately the skin friction attains the form (4.32). If the term 
}K were not present, this skin friction formula would be that appropriate 
to a steady free-stream velocity U,(1+-KT). The presence of the term }K 
means that the skin friction anticipates the free-stream velocity, since at 
time 7',, the skin friction has attained a state which, if the flow were exactly 
quasi-steady, would be appropriate at 7,,+-}. This anticipation arises 
because the pressure gradient which is accelerating the fluid has a greater 
effect on the more slowly moving fluid near the wall than on the fluid at 
greater distances from the wall. 
For steady flow with free-stream velocity U, the skin friction is 


T = pUy|vy| 








ee 
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but for large 7’ with free-stream velocity U,(1+AT7) the skin friction 
from (4.32) is ; ; ; 
Te. = ply Uw (1 +KT)+ tpl) V,.|K 
dU 
= plv,,|U(T)+4p\v,,| =- 
P Vy ( ) 4P Vu dT 


From (3.6) this may be written in the form 














. aU ma x 
T,, = plv,,|U(t)+ pd - = T,+ p04) —> (4.33) 
u PCy © (+ pojo dt a P10 dt 
where 8,, is the displacement thickness for a steady free-stream velocity 
ee 
4 
3 
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Fic. 3. Behaviour of skin friction in flow with a change of acceleration 


and 7, is the skin friction if the main stream were constant at U(t). It should 


be noted from (4.29) that as 7’ -> 00, 8,/v,,|/v tends to the steady form, 


namely unity; thus 5,5, . Consequently the displacement thickness 
approaches a quasi-steady state in which, although the fluid is accelerating, 
the steady displacement thickness applies. 

If the outer stream is at rest and the wall is moving then the force required 


to move the wall is given by the drag. The force required for a length / and 
unit breadth is thus 


5 ie us dU 
F = F,+ plbyo- = F,+pA 


dt” 


where A,, is the ‘displacement area’ (cf. Lighthill (1)). The term F, repre- 


(4.34) 
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sents the drag which would arise if the flow were steady at the instantaneous 
velocity U(t). We can thus call it the quasi-steady drag. The second term, 
however, is the product of a ‘virtual mass’, pA,,, with the acceleration, 
dU jdt. This ‘virtual mass’ arises because, as the wall is accelerating, fluid 
is accelerated as well through the action of viscosity (cf. Lighthill (1)). 


(6) Multiple accelerations 
From the results of section 4.3 (a) the velocity distribution can be easily 
derived for the case when the free-stream velocity is given a series of uniform 
accelerations. Consider, for example, the case in which the fluid accelerates 
from one velocity, U4, to another, U,(1+-K7%), between 7’ = Oand T = %, 
namely 


Uy (T <0) 
U(T)=({U(I+ KT) (O<T<N) 
U+KN) (1% <7) 
or f(T) = KT H(T)—K(T—T,)H(T—T,) 
so that {(7') may be written in the same form as (4.15) where here 
f*(T) = KT H(T). 
Precisely the same argument can be applied as in section 4.2 (6) to deter- 


mine the velocity distribution, skin friction, and displacement thickness; 
these are, respectively, 


u/U, = 1—e-7+a(n, T)—a(n, T—T), (4.35) 
where 
a(n, T) = 4K H(n)| sr—ar verte i 1") +-orfo( 7, de 1")! m 
4 afer nerfo( 7 “a m) a erfe( 4 m))|, (4.36) 


7,,|pUp|%-| = 1+A(T)—B(T—T), (4.37) 
where A(T) = 4K A(T) T+herf T!4+-2-*Te-T + T erf T'] (4.38) 


and 


8, |v |v = [L+7(7)—-y(T—T,) [14+ KT H(T)—K(T—T) H(T—T)), 


where (4.39) 
y(T) = K A(T)[a*Tte-7 +- 4a Te? —T*erfe T!+-T erf T'—}erf T*] 
for all finite 7’. (4.40) 


A graph of 7,, for K = 4 and 7, = } is given in Fig. 4. For T < N, 
A(T —T,) = 0 from (4.38), and (4.37) and (4.31) are identical so that, for 
T < }, Fig. 4represents the same skin friction as Fig. 3. As in section 4.3 (a) 
it can be shown that the skin friction decreases parabolically for small 
positive values of (77—7). For T > T, there is no pressure gradient so 
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that the flow will ultimately become steady; it can readily be shown that 
the ultimate velocity distribution and skin friction are given by 


u = U(1+KT)(1—e-) = U(o0)(1—e-") 
and 7,,/pU5|v,,| = 1+-KT, respectively. For K = 4, 7, = } these are 


u = 2U,(1—e-7) and 7,,/pU)|v,,| = 2. 
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Fic. 4. Behaviour of skin friction in flow with two changes of acceleration 


4.4. Decaying oscillatory velocity field 
Consider the case when the free-stream velocity is given by the real 
part of U(T) = Uj1—iA H( Tye}, 
where A, K, and Q are real constants. Consequently 
f(T) = —iA H(T)e“®* 107, (4.41) 


When K? < 1, f(7) satisfies the conditions of section 3 so that, from 
(3.12) the complex form of the velocity distribution is 


B= 1—e-1 iA HT) e107 4. 

U, 
4 (_ (na. 2) 

+iA H(T “ance TES | (K*-iN) exp} _ TT V3 
i (Te tat | e “7, (+75, ‘ da 
-} 4 
= 1—e-1~iA Hye feos a 
7 
0 


where a = (1—K?+iQ)! with the square root chosen so that re(a) > 0, 
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and where b = }». By use of equation (B, 2) of Appendix B, the complex 
velocity becomes 


. = ]—e-7—44 H( T )e “csr | — fer¥nehn Ks 4 
0 


x fe-na-K+400 epfo( 7 — (1K? 1Q)'TH) + 
l 47! 


+ erfe( +(1—K2+iQ)! "| (4.42) 


for T #0. Equation (4.42) has been derived on the assumption that 
kK? < 1 but it can readily be shown that u, as given by (4.42), satisfies 
(3.1) and the corresponding boundary conditions for all K? and Q, provided 
that 7 ~ 0; so that, exce,- at 7’ = 0, the velocity distribution corre- 
sponding to (4.41) is given by (4.42) for any real K and Q. The complex 
skin friction, 7,, may be obtained from (4.42) by differentiation with 
respect to » (ef. equation (3.15)); thus we obtain, for 7’ + 0, 
T(pUg\ty| = 1—}iAH(T)[et®* O78 + (77) -he-T + 
+e(K?-iOT (| — K®+ iQ) erf{T!(1—K?+iQ)'}]. (4.43) 
Consider firstly the particular case K2 = 1. Then the free-stream velocity 
is U(T) = Uf1+A H(T)e-? sin QT}, (4.44) 
and the complex skin friction is deduced from (4.43) as 
Tp pUy|t,| = 1—-iA H(T) e-7[ (2 T)-#+ OT + €'OT (iQ)! erf {((iQ7')4}]. 
But Jahnke and Emde (8) give the result, with a change of notation, 
erf {(ix)!} = (27)! C(x)—i S(z)], 
where C, S are Fresnel’s integrals, so that 
(¢Q)* erf {((tQ7)#} = (20) S(QT) +7 C(QT)). 
Consequently, when K* = 1 the skin friction is 
Trp pUg| tp 
= 143A H(T)e-7 [sin QT +(2Q)4S(QT) sin QT + C(QT) cos QT}}. 
(4.45) 
Consider now the case K2 4 1. When Q > |1—K?| the complex skin 
friction given by (4.43) is 
Tw/ PU! | 
= 1—}iA H(T)[(77')-te-P +e“ 107 + e117 (7) erf ((iQT')} | 
of which the real part is 
Tol PU) |X o| 
= 143A H(T)e-*"? [sin QT + (2Q)'{8(QT) sin QT +-C(QT) cos NT}}. 
(4.46) 


5092 .43 Y 
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This is identical with (4.45) when K*? = 1. Equation (4.46) may be written 
in the form 

7,,|pUy|%>| = 14+44e-™TH(T) B(Q, QT) sinf{QT+a(Q,QT)}, (4.47) 


where B(Q, QT) = [2QC%(QT)+{14 (20)! S(Q7T)}?}} 


(20)! C(Q7) | sii 
— “se . £ een ) 
(Q, OT) = tan 1+-(2Q)! S(Q7T) 





Thus (4.47) represents the skin friction exactly when A? = 1 and for 
QS |1—K?| when K? 4 1. Since « is positive, the skin friction has a 
phase lead over the free-stream velocity fluctuation. Note that, since 
O(x), S(x) > } as x > 0, for large Q7', (4.48) may be replaced by 


B(Q) = [Q+(20)'+1}) 


Q} (4.49) 
a(Q) = tan” 
and, if both Q7 and Q are large, (4.48) may be replaced by 
B= Qi, a = dr. (4.50) 


For any value of K®, the free-stream velocity is Uj{1+A QT H(T)] for 
very small 7’ so that the skin friction grows parabolically for small positive 
T (cf. section 4.3 (a)). As time increases the skin friction starts to fluctuate 
due to the fluctuating free-stream velocity and is given by (4.47). For 
sufficiently large values of the time, Q7' is large for any non-zero value 
of Q, and then B and « are given by (4.49) if K? = 1 or if QS |1—K?. 
If, furthermore, Q > 1, B and « are given by (4.50). In the latter case, 
the free-stream velocity and skin friction are given, for 7’ > 0, by 


U(T) = Uj1+ Ae-*’? sin QT] ) 

Tp pUp\Yyp| = 14+-}Ae-**T Yt sin(QT+}x) | 

But Stuart’s (2) result for high frequency may be written in the form 
U(T) = Uji+esin QT] ) 

Tp pUg|Y,| = 1+ 4eQ! sin(QT + 47) 

which formally corresponds to (4.51) if = Ae-**?, Thus the skin friction 
(4.51) behaves as though the main stream were a pure oscillation about a 
mean, with amplitude factor Ae-**?, By comparison with quasi-steady 
analogies obtained for other cases, we may refer to this case as quasi- 
harmonic. In other words this quasi-harmonic property follows from the 


fact that; for large Q, the logarithmic decrement of the second term of the 
skin friction (4.51), 


(4.51) 


(4.52) 


2. "2 
*log.¢ = 0434257, 


cy 
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is small. The skin friction given by (4.51) also formally corresponds to 
Lighthill’s (1) expression for the skin friction at large frequencies if 
« = Ae-*T, 

A graph of the (exact) skin friction for K = A = l and Q = 10is given 
in Fig. 5. The properties of the skin friction discussed above are exhibited 
in this graph. For large 27’, from equation (4.49), B = 3-93 and a = 0-605. 
It can be shown from the graph that B and a app‘owch these limiting 


un) “T 
Up 
cael 





© 1 it 1 5 











i 


I aii 
-0+4 o3 NOS 0506 07 08 09 1-0 





Fic. 5. Behaviour of skin friction in decaying oscillatory flow 


values more and more closely as 7’ approaches 1, so that the limiting values 
of B and « are attained for quite small values of 7. This is due to the large 
value of Q considered. The values of B and a which result from using the 
more approximate form (4.50) are B = 3-16 and « = 0-785, and these do 
not differ greatly from the values B = 3-93, a = 0-605. 
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APPENDIX A 
Suppose that the free-stream velocity ultimately becomes steady, that is, 


lim U(T) = U(oo) or lim f(T) = f(), 
T« Ta 
where U(co), f(o0) are finite constants. Then, from (3.12), 


2 


etn 


2 
lim (7) 1 e+ foo) — TF lim [ aty7—dyexp{—(a+7)} aa (A,1) 
"9 


T>@ 7 
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- 
Let x(A, 7) = A exp| — (a+). (A, 2) 


Then the integral in (A, 1) may be written in the form 


oe o A 
J S(T—A)y GA = f f(c0)y ad + 
0 


| {f(T—A)—f(ao)}x aA + | {f(T —A)—f(co)}x dA. 
0 Ay 

(A, 3) 
The suppositions that f(0) is finite and f(7’) is bounded in any finite interval 


imply that f(7') is bounded in —«o < T < o so that 





x a 
| {f(T'—A)—f(~)}x dA < A | A-te-A dA ¢s (A, 4) 
Xo Xo 
for Ay sufficiently large (A is a positive constant). Having chosen some A, to satisfy 
(A, 4), then, since | f(77—A)—f(o)| < ¢€ for T sufficiently large, 
Av As = 
| {f(T —A)—f(~)}x dA) < € | xdA <€| ydA= Ke (A, 5) 
| 0 0 


0 0 


for T sufficiently large, where K is a finite positive constant for 0 < » < «©. Thus 
from (A, 1), (A, 2), (A,3), (A, 4), and (A, 5) for 0 - 7<@ 


»>—4n P 
lim (~) = 1—e-+ f(«0)— 7 fl) | A-Pexp|—(A+ 
é 


T>« 0 


7? \\ 
1 an) - 


= l—e-"+f(0)—e-" f(a) = (l—e")[1+-f()] 
by equation (B, 4), Appendix B. 


Thus, from (3.12), lim uw = (l—e~”)U(o). (A, 6) 
) es 
Similarly, if lim U(T) = U(—o), where U(—oo) is finite then 
T>- « 
lim U(T) = (l—e)U(—o), (A, 7) 
T->— @ 


APPENDIX B 
Some integrals used in the text 


For re(a) > 0,6 > 0,2 > 0, 


Zz 


. 2 
I(x,a,b) = | exp{—(aa+5)}a4 dr 
0 
4 ‘b b 
— 7 abi p—4ab a be is =_ 
— Sa°” le 4 erfe(——ai#) erfc ([+ar')], (B, 1) 
r 


J(x,a,b) = 


. 


0 


° 2 
exp| _ (aA + 5) ha i dd 





4 ‘b b 
ie ee tab| esa ae eee fd 9 
= 35° ab/ ¢—4a erfe( ax ) erfe j tart) : (B, 2) 
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Particular cases of (B, 1) and (B, 2) are 


+ 
I(«,a,b) = —rm, (B, 3) 
+ 
J(0,a,b) = Femme, (B, 4) 
Furthermore, (B,1) is valid when b = 0 so that 
+ 
I(x,a,0) = ~erfazt (B, 5) 
al 
and I(~,a,0) = ae (B, 6) 
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SUMMARY 


Assuming that the solution is known to a three-dimensional elastic boundary- 
value problem for a particular Poisson’s ratio, equations are established which enable 
the general solution to be deduced from the known one. The known solution need 
not necessarily correspond to a real Poisson’s ratio; an infinite value, for example, 
reduces the problem to one in Potential Theory. As illustrations of the method three 
problems are discussed, one dealing with as 1 ‘nfinite body whose plane boundary 
is rigidly fixed, and two with cones loaded ut the vertex.. The solutions to the cone 
problems appear to be new. 


1. Introduction 


ALTHOUGH the mathematical theory of elasticity was firmly established 
by the end of the last century the small number of problems which have 
been completely solved testifies to the difficulty encountered in any 
practical applications. Even when a general class of boundary-value 
problems has been integrated the solution usually contains a set of in- 
tegrals which cannot be exactly evaluated unless the body has some 
elementary geometrical shape and is loaded in a simple or artificial 
manner. 

These difficulties are increased by the presence in the theory of the 
elastic moduli as parameters instead of as definite numerical quantities. 
Retention of the moduli as parameters ensures that the solution to a 
boundary-value problem is applicable to a general type of medium, but 
there are many instances where considerable simplifications in the working 
follow if numerical values for the moduli are inserted at the outset. The 
resulting solution is then applicable only to a single type of medium and 
so is strictly limited, but the question then arises whether the general 
solution to the same problem can be deduced from the known particular 
one. The answer evidently presents a new method of solving the elastic 
equations and in the present work the three-dimensional aspects of this 
topic are discussed. 

The two-dimensional case has been partly dealt with by Michell (1) 
(see also (6) p. 119), and Coker and Filon (2) who investigated solutions 
only when the boundary conditions were of the stress-type. The three- 

(Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 3, 1958) 
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dimensional problem was first considered by Westergaard (3) who sup- 
posed that, of the two elastic moduli, only Poisson's ratio was varied to 
give a simplification. We make the same assumption and follow Wester- 
gaard’s original idea and obtain the general solution from the particular 
one by the addition of a certain distribution of stresses and displacements. 
Westergaard’s approach made use of the so-called twinned gradient and 
was only partially successful. The presentation we give leads to a com- 
plete development, and suggests that there are several means available 
for determining the stresses and displacements to be added to the parti- 
cular solution. The one adopted here employs the Maxwell stress functions 
and includes Westergaard’s results as a special case. 


2. Preliminary theory 
It was shown by Maxwell (4) that for a medium not subject to body 
force the Cauchy equilibrium equations 


M4 = 0 (i,j =1,2,3) (1) 
ls 


referred to rectangular cartesian coordinates, have a solution 





, 8B, 80, _ 20, AA | _ ad OB 
= O22 * ayt? 4% Gxt " G22’ —* ay? * azt (2) 
0A eB e. eC 
Cc =, fog =—-—-——-, it detain 
” Gyez - Oz0x ”" exdy 


where the o’s are the cartesian components of the stress tensor, and 
., 8, C are arbitrary functions of x, y, z. The completeness of this solu- 
tion has been proved by Langhaar and Stippes (5). 
If the medium is now assumed homogeneous and isotropic then substi- 
tuting the stresses (2) into the Beltrami—Michell compatibility equations, 
viz. 





1 @o 1 @o 1 @o 
v? ee Se a a 
tian Metin Ost Thy oat 
1 @o l @o l @o 
Wie tet ee Deter ee We oe 
"wt ity aye?! Thy deen tity ony 


(3) 
in which o is the hydrostatic tension and v is Poisson’s ratio, shows that 


A, B, C satisfy the compatibility relations 


2A @B ec 
V24 — V2B = V20 = [-— 2—y), 4 
A B C (; st+a5 =) | v) (4) 








where V? denotes the Laplacian operator. 
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Also, from (2) the hydrostatic tension becomes 
2A &B a) 


c= (o,+0,+¢,) = V{A+B+0)—(—5 +3 22 


which, with the help of (4) may be expressed as 








a= (5*)p. (5) 
2—v 
2A &@B ec 
where D= tat + Gye + oat 


Adding the three Beltrami—Michell equations (3) containing the normal 


stresses gives Gs af (6) 


so that, by (4) and (5), 
V4A = VIB = VIC = 0; (7) 
that is, A, B, C are biharmonic functions. 


The displacements u = (u,v, w) are found by substituting the stresses 
(2) into the stress-strain relations 


» one v 
2MEey; = O45 — I+ 
where ¢,; are the cartesian components of the strain tensor, » is the shear 
modulus, and 6;; is the Kronecker delta. The resulting relations are 
simplified by means of (4) and then integrated to give the displacements. 
The additional functions introduced by the integration represent a rigid 


body movement and so can be omitted in the final form of the displace- 
ments, which become 


0835, 


eT a 
uu = £(4—B—0) 
é 
2uv = _ — 
Qu y B C—A) }. (8) 
Quw = 2 (C—~A—B) 
ae } 





It will be observed that neither of the expressions for the stresses or 
displacements in terms of Maxwell’s stress functions explicitly contain 
Poisson’s ratio. This fact is of importance in the present study. 


3. Effect of varying Poisson’s ratio 


Suppose that a given three-dimensional boundary-value problem with- 
out body forces admits of a simple solution when the Poisson’s ratio for 
the medium assumes some particular value v’, so that the Maxwell stress 
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functions A’, B’, C’ are known. Then in accordance with equations (2) 
and (8) the stresses and displacements may be expressed in terms of A’, 
B’, C’, and by equation (4) these functions themselves satisfy the com- 
patibility relations 

V2A’ = V2B’ = V°C’ = D’/(2—’) 
eA’ 8B’ aC’ . (9) 
(iar + et Ge) 

Now let Poisson’s ratio change to a new value v, but suppose that the 
shear modulus and prescribed boundary conditions remain as before. The 
elastic field is thereby altered, and usually the corresponding field equa- 
tions must be solved afresh. However, it is proposed to show that the 
altered field may be determined from that associated with the particular 
value v’ of Poisson’s ratio without any further explicit reference to the 
field equations. As such it appears to offer a new method for solving 
elastic problems in three dimensions. 

We will call the elastic field associated with the Poisson’s ratio v’ the 
initial state and designate with primes all quantities connected with it. 
The elastic field associated with the new Poisson’s ratio v and which we 
seek to determine, will be called the final state, and the corresponding 
quantities are left unprimed. This notation is adopted throughout. 

We assume that the final value of Poisson’s ratio is perfectly general 
and is only restricted by sufficient conditions for the uniqueness of the 
final state. For example, v must be within the range —1 < v < } in order 
that the strain-energy function be positive-definite, as required by Kirch- 
hoff’s uniqueness theorem. The theorem further requires the strains to be 
derivable from the displacements according to 


_ feu; , ou; 
«a= alee aa) 


where D = 


and to be related to stresses by means of the generalized Hooke’s law 
9 

oy = By) See st He (k = 1,2,3) (10) 
where the summation convention is used; the stresses, moreover, must 
satisfy the equilibrium equations (1). Besides these field equations, the 
elastic components must also satisfy the prescribed boundary conditions. 
Such a state, satisfying all the requirements of the uniqueness theorem, 
is called a possible elastic state and the associated Poisson’s ratio is said 
to be real. It should be noted that if the field equations are still satisfied 
for an unreal Poisson’s ratio the only consequence is that their solution 
is possibly no longer unique and therefore not physically significant. 
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The final state therefore forms a possible elastic state with a real 
Poisson’s ratio and the same boundary conditions as in the initial problem. 
A set of Maxwell stress functions A, B, C can thus be associated with the 
final state, and these unknown functions obey the compatibility relation 
(4) so that the final stresses and displacements can be derived from them 
by means of equations (2) and (8) respectively. 

Now consider the difference between the final and initial stress func- 


tions: Au hs, B= B-F, C= C—C. (11) 


It follows immediately from (7) that A”, B”, C” are biharmonic functions. 
By analogy with (2) we may formally define an equilibrium state of stress: 


»_ &B" 2c" » 80" A" » 4" FB’ 





Cc — CGC. —- —-_-— —_—— C.= —— 
z a8 I 7 ’ n.§ 4 ’ z a..5 
a2? * dy? . et " at dy? © oat re 
» (12) 
" 627A” ° 6? B” . e2C” 
Co = =_— mame 9 OC, == etna Co == ae 
- eyez r 20a — exoy 


and similarly by analogy with (8) we may formally define a set of displace- 
ments to be 





é 

9 ” an een eee Pes ” 

Quy 5, (A B"—C") 

” 0 ” ” ” 

Que” = —(B’—C"—A") ). (13) 
ey 

Su" = 2 (C0"—A"—B’) 
oz } 


The stresses given in (12) are not in general associated with the displace- 
ments given in (13) through stress-strain relations of the type (10), and 
indeed the precise connexion between them is by no means obvious. 
However, it is possible to interpret (12) and (13) in a very simple manner. 
For by the Principle of Superposition and equations (2) and (11) the stress 
distribution (12) forms the difference between the final and initial stress 
distributions and thus represents the increment in the stresses which 
occurs when Poisson’s ratio is changed in value. Likewise for the dis- 
placements (13); they may be regarded as forming the increments in the 
displacements of the medium due to a change in Poisson’s ratio. 

This interpretation provides the key to the present solution of the field 
equations, because once these increments have been determined the final 
state is essentially known. 

The most suitable way of calculating the additional state, as it may be 
called, is to derive equations for the additional stress functions A”, B”, C” 
in terms of the known components of the initial state, solve these equa- 
tions for the stress functions, and then use equations (12) and (13) to 





| 
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a al aaah 5 
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obtain the additional stresses and displacements. The equations for A’, 
B", C” are established below, but before doing so we briefly mention 
another approach which although not so convenient in the solution of 
actual problems is nevertheless useful in the deduction of certain theorems 
concerning the additional state.t This second method consists in deriving 
equations for the additional stresses and displacements directly by form- 


ing the difference between the appropriate field equations of the final and 
initial states. For instance, since 


O;; = O3;—G;; 
u” = u—u’ 


we may take the difference of the respective stress-strain relations (10) 
and obtain 


, 
0 nae 


” v ” ’ ” 
oi; = - taet ji 8+ 2H (14) 





= (1—2v)(1—2v 


which relates the additional stresses and strains. Again, those equations 
which are derived from the Beltrami—Michell and Navier compatibility 
relations may be identified with those of the classical field theory, but 
include a body force term whose magnitude depends upon the initial 
hydrostatic tension and the initial and final Poisson’s ratios. The solution 
to these equations presents no essentially new approach and for this 
reason we omit a more complete discussion. 

Any calculation of the additional state is simplified by the boundary 
conditions for this state, since these are homogeneous. That is, on those 
parts of the surface where the tractions are prescribed the additional state 
has zero tractions; on those parts of the surface where the displacements 
are prescribed the additional state has zero displacements; and on those 
parts where the boundary conditions are of the ‘mixed’ type, the additional 
state has the corresponding ‘mixed’ components zero. When dealing with 
the stress-functions A”, B”, C” it must be remembered that the boundary 
conditions involve their first and second derivatives, and not the functions 
themselves, so that it is appropriate combinations of these derivatives 
which become zero on the boundary. 

Equations for A”, B”, C” may be derived as follows. Subtracting (9) 
from (4) and using (11) gives 


V2A” = V2B" = V2C" (15) 


D D’ 


20" — 7 J 
and also Vv2C @—)~ @—r) 





(16) 


+ A full description of these theorems is beyond the scope of the present study. 
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Substituting 
, ‘ as 027A” @%B" aC" 
D=D+D, D'= ( a a 


in (16) we obtain (2—v)V20"—D" = (= \. 





, 
2—y 


which on using (5) can be rewritten in its final form 
eA” 3B” 3C" v—v’ 
a pelle ll eta FO states 17 
a—wwor— (Fa + oe + Se) = (ie) (17) 
Equations (15) and (17) together with the boundary conditions enable 


the biharmonic functions A”, B”, C” to be determined. This is illustrated 
later in sections 8 and 9. 


4. Axially symmetric states 

For many practical applications it is useful to have the results of section 3 
reformulated in the special case of axial symmetry. We use cylindrical 
polar coordinates (r,6,z) with z as the axis of symmetry. Langhaar and 
Stippes (2) have shown that the complete representation of the stresses 
now requires just two biharmonic functions F, H, depending only on r 
namely, 

a _ @#F 10H al OF | eH i. _ OF Ler 


—y 5 


~ 
> a 








tO" ¢ Or’ éz? " Or? er?" ¢ Or 
(18) 
Cc. = — “ 0,9 = = 0 
 Oré2’ tts. 28 
The Beltrami-Michell compatibility equationst 
2 1 @o 
v2 its ee ome: auaenios ates Se @ 
Or Far 90) + ope 
2 1 léo 
V2 ys ome ht atmimny «scum S08 0 
Ot Fa (80) + Te oe 
a (19) 
V%.4+.—__.. __ = @ 
Tf») at 
1 1 ea 
V°0,,—=°%: +-——. —— = 0 
re 8° Dy) roe 
after substitution of (18) lead to 
or — v2 ] a? . 
V2F = V?H = —__ © (F—H). (20) 


(v—1) éz? 


t The compatibility equations in the axial-symmetric case are derived in (6), chap. 
xiii. 
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From (18) we have for the hydrostatic tension 


e2 
o= VF +H) + —(F—H), 


pe 
so that with (20), c= ot ne F—R). (21) 
v—1] 2? 
The displacements (u,, ug, w) associated with the stresses (19) may now be 
calculated as in section 2. On integrating the stress-strain relations after 
these have been reduced by means of (20) the displacements to within a 
rigid body movement become 
2uu. = _ eH 
iid ér 
ug = 0 ; (22) 
. e 
2uw = —(H—2F) 
z 


To calculate the final state from the known initial one we proceed 
exactly as before and assume that (F’, H’), (F, H) are the initial and final 
sets of the Maxwell stress functions. Then F’, H’ satisfy 


Vir’ = Vin’ = |. 
(v —1) 


with similar relations for F, H. The stresses and displacements corre- 
sponding to both states may be derived from equations of type (18) and 
(22). By taking 


e , , 
aga —H’) (23) 


F" = F—F’, H” = H—H' 


we see that the stresses and displacements defined by equations (18) and 
(22) in terms of F”, H” represent the increments in the corresponding 
components due to a variation of Poisson’s ratio, and hence are the stresses 
and displacements of the additional state. Subtracting (23) from (20) 


aia vee" — V2H" (24) 


and o—1yver’— &(Fr_H") = ()e" (25) 
Oz l+y 
which are analogous to equations (15) and (17) of section 3. 

The boundary conditions for the additional state are again homogeneous 
and similar conclusions follow as in section 3. Using these boundary 
conditions and (24) and (25) we may determine F”, H” and thus the 
additional stresses and displacements, and then complete the problem as 
previously described. 
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5. Connexion with Westergaard’s method 
If in the general case of section 3 we put A” = B” = 0, and in the axial- 


symmetric case we put F” = 0, both C” and H" become harmonic func- 
tions. Also (16) and (25) reduce to 





<p (=o as) 
ez? oz" i+v 
Moreover, from (16) we have 
D D’ 
(2—v) (2—»’)’ 
which on using (5) becomes 
eas Sih setae (27) 
l+y l+v 
and so o” =o—o' = (=) (28) 
l+yv 


Similar results hold for the axially symmetric state. 

Equations (26), (27), (28) have been previously deduced by Wester- 
gaard (3) by an allied method which makes use of the so-called twinned 
gradient. In this method Westergaard restricts the additional displace- 
ment field to that which can be expressed as the twinned gradient of some 
function ¢ to be determined later. That is, 


eet ,_o ,8 o e 
2uu” = (15-15 +k 5 }¢, (29) 
where i. j, k are the unit cartesian vectors, and the operator on the right 
side is the twinned gradient. This form of u” is clearly axially symmetric 
about the z-axis. The additional strains are derived from the displace- 
ments in the usual way and then the stresses are obtained by means of 
the stress-strain relations (14), although Westergaard deduced these in a 
less straightforward manner than that given here. The additional stress 
distribution so obtained is also axially symmetric with the shear com- 
ponents o,,, o,, both zero. In order to establish equations which can 
easily be solved for ¢, Westergaard had further to restrict this distribution 
to be planar by taking o, zero. Then ¢ becomes a harmonic function 
identical to C” and H”’, and satisfies equation (26) which can now be 
solved. The method is applied to the problem of a force at a point of the 
plane bounding a semi-infinite body and also to the problem of a thick 
disk rotating about its axis of revolution, but in both these cases it so 


happens that the additional state is a planar one, otherwise the method 
would be inapplicable. 
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Westergaard’s approach neglects any mention of boundary conditions 
for the additional state but inspection shows that he assumed these to be 
of the ‘stress-boundary’ type. A further requirement, therefore, is that 
the additional surface tractions must be zero. This is automatically 
satisfied in the case where the body is bounded by planes, the additional 
state of stress occurring parallel to these planes. 

The restrictive character of the method is caused by adopting the form 
(29) for the additional displacement instead of calculating it as we have 
done. Indeed, although Westergaard includes in his presentation an 
indication for a generalization,+ this artificial step obscures the under- 
lying principles and it is only when we regard the additional state as 
the difference of two possible elastic ones that an extension becomes 
possible. 


6. Choice of initial state 


In choosing the initial state we are not prevented from considering any 
value whatsoever for the initial Poisson’s ratio, since although the field 
equations remain physically significant only for real values of v, they are 
mathematically consistent for ali values and hence yield a solution no 
matter what value of Poisson’s ratio is taken. Thus, any boundary-value 
problem can be solved with an unreal v, and moreover the same methods 
can be adopted as for the ‘real’ theory, except that now the ‘stresses’, 
‘strains’, and ‘displacements’ are purely mathematical entities and the 
solution is possibly no longer unique. The derivation of equations (15) 
and (17) depends only on the mathematical validity of the field equations 
and so is true for all values of v and v’. 

Hence we may solve the given boundary-value problem for an unreal v 
and use our method to convert the solution to one corresponding to a real 
value. The final solution is unique since all the conditions of Kirchhoff’s 
uniqueness theorem are met. 

The advantage of extending the choice of initial Poisson’s ratio to unreal 
values lies in being able to select any value which offers a simplification. 
Apart from values which simplify particular problems, there are also 
others which when substituted into the field equations effect a significant 
reduction to the theory as a whole. Examples of these values are 0, 3, 
+1, «©. The real cases 0, } have been dealt with by Sokolnikoff (7); we 
consider in the next section the unreal case vy = 0, and then use the results 
in the solution to a particular problem. 


+ The present extension of Westergaard’s original idea was first suggested to me by 
Professor R. Hill. 
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7. Field equations for infinite Poisson’s ratio 


When Poisson’s ratio is infinite and provided the shear modulus is held 
finite, Young’s modulus is also infinite, and Lamé 3 constants A, u are 
connected by the relation 

A+p = 0. (30) 


The stress-strain relations (10) become 


Oj; = p(2€,;—€ 44 5;;) (31) 


and hence by the addition of the normal components of (31) we have 


Cc => — BE pp 
From (1) and (31) with the strains expressed in terms of the displacements, 
V-u = 0. 


That is, the displacements are harmonic functions; and therefore so also are 
the stresses and strains. Moreover, the equilibrium equations are still 
identically satisfied by (2) and similar arguments which led to (4) now 
show that the Maxwell stress functions, instead of being biharmonic, are 
harmonic. 

The solution to Kelvin’s problem of a point force inside an extended 
body becomes particularly simple. Suppose the force is of magnitude P 
and acts at the origin of coordinates along the positive direction of the 
z-axis. The specification of a point force requires the stresses to be O( R-*) 
near the origin, where R is the distance from the origin; also, for an infinite 
medium the stresses must vanish as O( R-*) at infinity in order to comply 
with Kirchhoff’s uniqueness theorem. These conditions apply for real 
values of Poisson’s ratio but we assign them to the present problem so 
that the stresses of any derived elastic state containing a point force have 
the correct order of magnitude. 


To describe completely the stresses we need only two stress-functions 
as there is symmetry about the z-axis. These functions must be harmonic 
and give rise to stresses of O(R-*) both at the origin and at infinity. 
Therefore, let ut F = Blog(R+2), 


where B is a constant to be determined later. 
In accordance with (18) and (22) the stresses and displacements are 
== — va ’ ; ... = “ = Ni 
0, q» 9G Be o, Bas Ors Br (32) 
u, = Ug —— 0, 2uw = —2B) R 


To find B we observe that the total vertical pressure on any plane 





ae OO 
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z = a positive constant must be }P. Hence, for z > 0, 


\P = —2nBz i 
and so B= - 3 
4n 


Substituting B in (32) completes the solution. 


8. Internal force in semi-infinite medium with plane rigidly- 
constrained boundary 

We now use the solution to Kelvin’s problem to solve the problem of 
a force inside a semi-infinite medium whose plane boundary is held rigidly 
fixed; that is, all components of the displacement are zero on the boundary. 
For an infinite Poisson’s ratio the derivation from Kelvin’s solution corre- 
sponds to a problem in Potential Theory, and can therefore be easily 
accomplished. We thereafter use the results of section 4 to obtain from 
this particular state that associated with a general value of Poisson’s 
ratio. 

Let the medium occupy the space z < 0 and take the force to be of 
magnitude P acting in the positive direction of the z-axis at the point 
(0,0,—c), Due to this force the solution contains a singularity of the 
Kelvin type found in (32), and further this can be the only singularity 
since there is only one force acting. Therefore if S; represents the distri- 
bution (32) with asuitable change of axes, we may assume that S, provides 
the singular part of the final distribution and so when Poisson’s ratio is 
infinite the problem is reduced to finding a distribution S, regular and 
harmonic in z < 0 such that the displacements belonging to 

S =8,+8_ 
vanish on the boundary z = 0. 

Now in (32) the w-displacoment is clearly the singular part of the first 
Green’s function for the half-space occupied by the body. This function is 
harmonic, vanishes over the boundary and is, moreover, known for the 
region under consideration. Hence it will serve for the w-displacement 
of S. The u,-displacement of S, is identically zero throughout the body 
so that the corresponding displacement of S, must vanish on z = 0. Since 
this component is regular we have from a well-known theorem in Potential 
Theory that it is identically zero throughout the half-space. Thus we 
may write for the displacement of S 


u, = 0, 2pw’ = wla- ze) 
1 


5092 .43 Z 


(33) 
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where R? = 2?+y?+(z+¢)?, R? = 2?+y?+(z—c)?. 


The second displacement component, ug, is obviously zero since there is 
axial symmetry about the z-axis. The strains are 


eee 8 é = al I 
“* 4mpez\R_ Ry’ “Sap or\R R,) J, (34) 
€-=@= 0 


and from (31) the stresses are 


— ae. 2 aa. a 
Or ~~ Ge z\R ral 0 an aa\R_ R, 


Pee ck. | Tee” Pi 2) | ee 

° 40 é2z\R- R,)’ " 4nér\R R, 

Equations (33), (34), (35) supply the solution to the problem when 
Poisson’s ratio is infinite. 

The solution for any other value may now be derived from the theory 
described in section 4, with the initial state corresponding to one with an 
infinite Poisson’s ratio. Thus, we must find two biharmonic functions 
F", H" which satisfy v2F" — V2H" (36) 


(35) 


a2 
and (v—1)V2F"—<_(F"—H") =o’, (37) 
C2 


where o’ is given in (35). Furthermore, since the displacements are every- 
where prescribed on the surface, the boundary conditions for the addi- 

tional state are that the displacements vanish everywhere on z = 0. That 
i u.=w’=0 onz=0, 

0H” é 

or by (22), — = —(H"—2F")=0 onz=0. (38) 
or ez 


To calculate F”, H” we start by integrating (36), giving 
H” = F"-+a, (39) 


where « is a regular harmonic function in z < 0, and then substitute for 
H”" from (39) in (37). This gives 


w. ty» 2x 
VF = {0G} /[o—». 


which on using (35) for o’ becomes 


Vir" — |-z é i ae, _ a (v—1), 
4ndz\R R,} 22? 
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and so, by integrating, 
1 zte 2z Zz Ox 
F’ = Sone Abe sae siping pcos 
staal R z 2(1—v) a t Bs 


where § is also a regular harmonic function in z < 0. Returning to the 
boundary conditions (38), we put 





(40) 


H"=0 onz=0 (41) 
and thus obtain by (39) 
F”’=-—a onz=0 (42) 
and by (38), il — ~ on z = 0. (43) 
éz ez 


Then since « is regular we have from (42) and (40) 


Pe 1 " nt 


and from (43) and (40) 


; i—S\ea Of) 
l—v})é@z = 8a(1—v) ez R,}  @ 


which, on integrating once with respect to z, becomes 


1 /1—2v P c 
3° i) ” -saale)t B. (*) 


No additive arbitrary function of x, y is required since «a, 8, 1/R, vanish 
as z—> —0O. 
From (44) and (45) we-find 


a (i) 
~ 3n(3—4y)\R,)’ 


P c 
dus saaagacale 


and these expressions when substituted into (40) and (39) lead to the 
following final forms of F” and H”: 


PF = PFI /[teeu—n+ Pel ee -25(F -)| [tas 4v)(1—v) ], 
HH’ = Pleto(a—z ]/texa—n)]—Pee Ep) [[4n3—40)0—»)} 


(46) 
Calculating the additional state from (46), (18), and (22) and adding this 
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to the stresses and displacements of (35) and (33) respectively gives the 
final solution to the problem for any value of Poisson’s ratio: 


2d Ee a om a 
4 slp) — mam RS Rs Re RS 
essen [3 a) 2a (i 
82(3—4v)(l1—v)| e2*\R, ér?éz\R,) | 


= Pr of} 1 ‘a Pp z—c ee |A 
_ -Zale- zi) '87(1—v)| R RF RS 


3Pc fe=e-20 1 ' 


o> — 





~ 8x(3—: 4v)(1—v) R} R 


ros) l Pp ee zte  32%(z—c) | 
— ~—- =. fae: + 
‘Sale R, 87(1—v)| 


Re R32 OR Rs 
rye x) 
' 82(3—4y)(1—v)| éz3 \R,J ez? R,) |’ 


Pe 
_Paji 1 P 11) _3(+e)*r | 3rz(z—c) 
n= aR 5) aaa’ Re at 


1 


> 
7 














Pe 3r(z—c) 6rz __ 42r2(z—c)? 
+Sasagaca RS ‘RS BR "|; 
— -c(E)ltecok ail tea ain) — 
tia |" + cane 
= a5" e| tae ey 


om E(t 1)4 FP _[teteh_ toh, 
27\R R,} ° 8x(1—v)| RB R' 
Pe _[3(z—c)? 1 
+m Rl 

The solution was first obtained by Rongved (8) who used a slight modi- 
fication of the method adopted by Mindlin for the allied problem of a force 
inside a semi-infinite body whose plane boundary is traction-free. Both 
approaches are closely related to Potential Theory and therefore require 
a knowledge of the Green’s functions for the region occupied by the body. 
They cannot be generally applied and are successful only in those prob- 
lems where the Green’s functions are known. Our method is of course not 
bound by this restriction since it is only for an infinite Poisson’s ratio 
that Potential Theory is explicitly introduced. However, it is preferable 
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to use Mindlin’s approach since this gives immediately the final solution 
for any value of v and does not involve the two subsidiary problems which 
are an essential part of our method. 


9. Two further problems 

To illustrate another way of determining the additional stress functions 
we solve the two problems of a solid bounded by two right circular coaxial 
cones loaded by a force acting at the common vertex in the directions 
(i) along the common axis; and (ii) perpendicular to the common axis. 
We suppose those parts of the solid at great distances are held fixed. 

By considering the form of the initial hydrostatic tension it is possible 
to express the additional stress functions as linear combinations of certain 
harmonic and biharmonic functions and thereby reduce equation (17) to 
a relation between the constant coefficients of the linear combinations. 
These constants are adjusted so that the reduced equation (17) and the 
boundary conditions for the additional state are satisfied. 

Westergaard, amongst others, has mentioned (3) that when Poisson’s 
ratio is } the stress distribution in Kelvin’s problem is purely radial and 
hence the stresses, strains, and displacements in this case are the same as 
those for variously shaped cones loaded by a force at the vertex, and also 
for an infinite wedge loaded at a point of its edge. Included amongst these 
general bodies are those considered in (i) and (ii) so that their solution is 
known for Poisson’s ratio 4 and we therefore take this to be the initial 
state. 


(i) Solid bounded by two coaxial cones loaded at the common vertex by a 
force acting along the axis 

Let the common vertex be the origin of coordinates and draw the 
positive z-axis coincident with the common axis of the cones. Take the 
semi-vertical angles of the outer and inner cones to be a and f respectively, 
and suppose that P is the magnitude of the force which acts at O in the 
positive direction of the z-axis (see Fig. 1). 

It is convenient to use the spherical polar (R, 6,4) components of the 
stress tensor and displacement vector and to express these components 
in terms of the coordinates (r, R,z) where 


r= Rsind, z= Reos¢. 
When Poisson’s ratio is } the solution may be derived from Kelvin’s 
problem and takes the form 


‘ . 3P (z 
on = Co => -3(R) (47) 
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with the other stress components vanishing; and 


r 


Zz P 
utp = -— —, =0, uu,=— ———, 48 
HOR sa RO wig ona R2 i) 


where a = (cos* B—cos* a). 


It may be easily verified that the stress-distribution (47) leaves the sur- 
faces of the cones ¢ = a, ¢ = f free from tractions and exactly counter- 
balances the force at the vertex. Any stress distribution which is added 








to allow for a variation in Poisson’s ratio must therefore leave the surfaces 
stress-free and also give rise to no further resultant force at the vertex. 


Clearly there is symmetry about the z-axis, so we may use the results 
of section 4. With 


3P oz 
2na R?’ 


equations (24) and (25) for the additional stress functions F”, H” become 


v’ = t, o’ = 


2 ~ 
(—1)v8F"— SF" H") = iow. : (49) 


2na RS’ 
ViF" = V8H’. (50) 
In order to reduce (49) to a relation between constants and also to satisfy 
(50) it is reasonable to try the following forms for F” and H”’: 


FP" =x k= +nlog(R+2)-+mlog( R—z) 
: (51) 
H’ = k= +hlog(R+2)+glog(R—z) 


where g, h, k, m, n are constants. Substitution of (51) into (49) gives 


2y—1)k-+ (hg) +p = (1-29) (52) 
with p = (m—n). 
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The cylindrical polar components of the additional stresses may be 
obtained from (51) and (18) and then transforming by the tensorial law 
we find that the stresses in spherical polar coordinates are 





























. ‘ or h g sin’d 
a= ee Freer: eal | 
in cos¢ _ h(cos? ¢-+-cos ¢—1) (cos? 6—cos d—1) 
4-0" R R*(1+-cos ¢) +9 R*(1—cos¢) 
ne oat + h g cos? d . (53) 
09 = (P— kh) +| sos t Toosa| 
tt h a dich 
ong = (P—F) +lizews =r 
os = O49 = 0 
Similarly, the displacements are 
. ost _ h g sin? ¢ ) 
te = Eee fee eal R 
us = 0 (54) 
a ES = h g sin ¢ 
anal la + b)-+0084 (+ salle ) 





The condition that the additional stress distribution (53) gives rise to 
no extra resultant force at the vertex is that the tractions in the direction 
of the axis of the cone taken across any spherical surface, centre the 
vertex, must be zero. That is 


| (ong sin? ¢—o', sin ¢ cos 4) dd = 0 
8 








or from (53), k{cos? «+-cos a cos B-+-cos*B—1]+p = 0. (55) 
The sides of the body remain stress-free, so that 
0% = Ong = 0 for ¢=a,f8 
and so from (53) 
‘ heos « g COS « 
(p—k)+ os - 
(1-+cosa) (1—cos«) (56) 
hcosB gcosB 7” 
(P—)+ (eos) * (1—cos8) 
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Solving equations (52), (55), (56) for the four unknown constants g, h, k, p 
gives 
(1—2v)P L(1—cos «a)(1—cos 8) \ 

















= 4ra M 

ee (1—2v)P L(1+ cos «)(1+-cos B) 
4a M 

ee (1—2v)P cosacosB 
dares M (57) 
_ __ (1—2v)P (L—1)cos «cos B 

il 4ra M 

where 
L = (cos* «+ cos « cos B-+-cos? B) 
M = (cos? «+ 2v cos « cos B+-cos? 8) } 


Substituting these values in (53) and (54) and adding the resulting system 
of stresses and displacements to the initial state, (47) and (48) give the 
final form of the distribution for any value of Poisson’s ratio 

P (83M+(1—2v)(L—3)cosacosB] z 


= — 


Qn LMN Rr 








(1—2v)P[ (1+ cos «)(1+ cos B) Fen en 
4n MN (R+2z) = (R—z) R3’ 
(1—2v)P 212R— R?) 


(2 
% = TN | (1-+008.) 1+-cos B) (R42) ~ 


—(1~e08.)(1—co8p) “EFF _ 2.08088 75, 


_ (1—2v)P 22 
°% = | (1+208y x)(1-+-cos f) ———. R(Ri2)~ 


—(1—cos «)(1—cos f) RS" 2 cos a cos 8 ah 


org = — = yr Ty | (1-084) (1-++-cos B) 





RTA 


—(1—c08)(1—e08 8) Far-R—> BR 2008 208 B Z| 
Cre = O46 => 0, 
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P 


= “Guy |——4 a)(1-+cosB) (1—cos | 








4nMN (R+z) (R—z) 
Ug = 0, 
2uug = to 2{L+-cos «cos B(3L—1)}—~. + 
$ 4n LMN RS 


v2 


{ Ne ax rz 
eRe here) |) ihe cos a)(1 008) aR | 
with ZL and M given by (57), and 


N = (cos B—cosa). 


The solution to the solid cone loaded by a force at the vertex along 
the axis may be deduced by taking 8 = 0, The problem has been pre- 
viously treated by Michell (9),¢ who combined various simple solutions, 
all containing a singularity in the stresses of O( R-*) at the origin, so as 
to satisfy the boundary conditions. 

By putting « = 47 in the solution to the solid cone we can obtain that 
for a force acting perpendicular to the plane boundary of a semi-infinite 
body at a point of the boundary. The solution was first derived by 
Boussinesq,t and later used by Westergaard to illustrate his method of the 
twinned gradient. Westergaard took as initial state Kelvin’s solution 
when Poisson’s ratio is 4, and indeed, had we directly applied our method 
to the same problem the various steps of the derivation would have been 
precisely similar, showing that the twinned gradient approach is a parti- 
cular case of our own. (For details of Westergaard’s solution, see (3) 
p. 139.) 


(ii) Solid bounded by two coaxial cones loaded by a force at the common 
vertex perpendicular to the axis 

We use the same general approach as in (i) with the exception that as 
there is no axial symmetry we must employ the results of section 3 instead 
of section 4. 

The common vertex of the cones is again taken as the origin of coordi- 
nates and the axis of the cones as the positive z-axis. The x-axis is drawn 
along the line of action of the applied force which is of magnitude P, and 
a and £§ are the semi-vertical angles of the outer and inner cones respec- 
tively. (See Fig. 2.) 

The initial state is chosen as the solution to Kelvin’s problem with v 


+ See also Love (10), 201. 
} Ibid. 189. 
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as 4, since then the stress-distribution is entirely radial; in cartesian 
coordinates this is 


, SP i{# ee: xy? of = —2P (x2 
a -=(m) y— wrab\ RS)’ -  rab\ RS 


» (58) 
» _ __ 38P [ayz » _ __ 3P (xz ay _ 3P (x*y 
Oey rab\ RS)’ 9 ab RS)’ rab VRS 








1x 
Fic. 2 


where a = 3—(cos*?«+cos «cos 8+cos? 8), 


b = (cos B—cos«). 


The initial displacement field is 


' P a*\ 1 : P [xy . P [xz 
2uu’ = — (1+ —,)— Que’ = — | Quw’ = (~~). (59 
“a ail +a - ~aaR a =a) (59) 


The stresses (58) leave the surfaces of the body ¢ = a, ¢ = B free from 
tractions and counterbalance the applied force. The additional stresses 
must therefore leave the surfaces stress-free and give rise to zero resultant 
force at the vertex, and so we adjust accordingly the constants of the 
linear combinations for A”, B”, C”. On substituting 





ied PP ..  S 
s ~ rab R8 

in equation (17), the relations for the additional stress functions A”, B’, 

C” become 





eA” eB" oF) -hne) 
3]? 


2—v)V20"— me tenn | ba 
(3-0 — Ft et ae) oe (00) 


V8A* == V2B" = V8C", (61) 
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Equation (61) is identically satisfied and equation (60) reduces to a relation 
between constants if A”, B”, C” have the form 





A’ =k ier Riz :)+ nl = *_)+alog(R+-2)+glog( Ra) 


B’ = k— nt? pes) +m( Rg) tmlog(R+-2)+-flog( R—2) (62) 





C” = k= Rt Macs 5) +a go) +m log R+2)+ flog R—2) 


where d, f, g, h, k, m, n, p, q are constants. 

With these expressions the cartesian components of the additional 
stress-tensor can be obtained from (2) and by transforming into spherical 
polar components it can be shown that the condition for stress-free 
boundaries requires the constants to satisfy the following equations: 























g=q=9%, (63) 
ei d cos* x 2p cos « 2ncosa __,) 
(k+p+n4 + Teo at (l—cosa)? (1+ cosa) fi-ood"" 

i cl a eer d cos? B 2p cos B 2ncosB _ 

(k PT N+) + eos Ap (1—cosB)? (1+ cos) ‘ (1—cosB) 

’ h cos «(1—cos? a) d cos a(1—cos? x) _ 

(k+p+n-+t)cos « — ia eae m 
p(2cos*a—1) , n(2cos*a—1) ; 

~~ (1-+e0s a) (l—cosa) 

’ hos B(1—cos*B) _dcosB(1—cos*B) _ 
Ctet+e+hemr— (1+-cos B)? ve (1—cos B)? 
__p(2.cos? B—1) n(2 cos B—1) _0 
(1+-cos B) (l—cosp)  ~—iyW 
where t = (m—f). 
Substitution of (62) and (63) in (60) gives 
2k(1—v)-+(h+d)—(p-+n)+t = (2v—1) =. (65) 


For the additional stress distribution to have no resultant force at the 
vertex, the tractions in the z-direction over any plane z = a positive 
constant must be zero. Thus 

27 ctana 
or drdé = 0 
0 ztang 


; » eB" 
or since o 


oo os then (a—2)k—p—n—2t = 0. 
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Solving the six equations (64), (65), (66) for the six unknown constants 
d,h, k, n, p, t gives: 

















d = _(1— 20)’ (2p cos B+ Bos.) | 
Abe 
2 R2 
h= n-s92- A?’ B*(C cos B+-D cos «) 
4n be 
: he 
ee), Ye 
ee a abe 
; (67) 
P CDE cosacosB 
—— aw 9 nol 
oe "om be 
P AB(2—E)cos «cos 8 
——- —— ae 
a 2a be 
— ~tia)* (aE cos «cos B—s) 
7 abe ) 
where 
a = 3—(cos?a+cos a cos B+-cos? f) \ 


b = (cos B—cos «) 
1 = cos? a-+-cos* 8—cos B(1-+-cos a cos f) 
8 = cos a cos B{1+-cos a cos B(cos « cos B—2)} (68) 
e = (1l—2v)s—al 
A = (1+ cosa), C = (1—cos«) 
B = (1+ cos 8), D = (1—cosf) 
E = (1—cos «cos 8) ) 








The final stresses «rd displacements for any value of Poisson’s ratio 
are then found by adding the additional stresses and displacements found 
from (12), (13), (62), (63), and (67) to the corresponding initial components 
in (58) and (59). This leads to 


P [3la* | (1—2v){ - 
Ee =. a COS Boat 


22 gp inte " 
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with a, 6, l, s, e, A, B, C, D, E given by (68). 
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We may again obtain the stresses and displacements for the solid cone 


by putting 8 = 0. The results have been previously presented by Michell 


(9) who used the same approach as in the cone loaded by a force along the 
axis. 


Further, by putting « = 47 in the solid cone solution we can obtain 


the solution to Cerruti’s problem of a semi-infinite body loaded by a 


ta 


ngential force at a point of its bounding plane. Westergaard illustrated 


his twinned gradient method by this problem. 
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SUMMARY 


A mathematical model is set up which describes the dynamical behaviour of a 
bowl which rolls on a deformable plane surface. Approximate solutions of the 
equations of motion are obtained, adequate for prediction of a bowl’s behaviour 
under normal playing conditions. These give the equation of the path and show that 
the angle of deflexion from a straight line is independent of the velocity of projection 
provided the bowling green conditions rem.sin unchanged. 

The intrinsic equation of the bowl path is found to be of the form sa 1—e-™, 
where the parameter p and the constant of proportionality depend on the type of 
bowl and the green conditions. 

The results are compared with observations on bowling greens and test-tables. 


1. Introduction 

WHILE it is well known that the curved trajectory of a bowl results from 
the superposition on its forward motion of a precession caused by the bias 
of the bowl, so far no complete dynamical investigation appears to have 
been given (Walker (1)). Since observation shows that there is no slipping 
between bowl and rolling-surface the assumption of retardation by a 
frictional force proportional to the normal reaction is precluded. The 
analysis here presented assumes that the effect of the rolling-surface in 
retarding the forward motion can be represented by a couple K, of constant 
magnitude K, acting on the bowl, the axis of K being parallel to the ground 
and normal to the direction of motion. Comparison with observation has 
shown that the resulting mathematical model is adequate for a wide range 
of initial and rolling-surface conditions. In particular there is close agree- 
ment between theoretical and observed behaviour in the game of bowls. 
Some details of the comparison are given in section 10. 


2. Frames of reference 
Fig. 1 shows the axial frame Oxyz, fixed on the horizontal rolling-surface. 
O is the initial point of contact of the bowl with the surface, Oz is the initial 
direction of motion and Qz is vertically upwards. The distance OP 
travelled along the path in time t is s, the range of the bowl is R, and the 
total angle of deflexion A. The coordinates of the centre of mass @ of the 
bowl at time t are (z, y, 2). 
A vertical section of the bow] is shown in Fig. 2. The moving axes Gn, 
[Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 3, 1958) 
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SUMMARY 


A mathematical model is set up which describes the dynamical behaviour of a 
bowl which rolls on a deformable plane surface. Approximate solutions of the 
equations of motion are obtained, adequate for prediction of a bowl’s behaviour 
under normal playing conditions. These give the equation of the path and show that 
the angle of defiexion from a straight line is independent of the velocity of projection 
provided the bowling green conditions remain unchanged. 

The intrinsic equation of the bowl path is found to be of the form s « 1—e~?¥, 
where the parameter p and the constant of proportionality depend on the type of 
bowl and the green conditions. 

The results are compared with observations on bowling greens and test-tables. 


1. Introduction 

Waite it is well known that the curved trajectory of a bowl results from 
the superposition on its forward motion of a precession caused by the bias 
of the bowl, so far no complete dynamical investigation appears to have 
been given (Walker (1)). Since observation shows that there is no slipping 
between bowl and rolling-surface the assumption of retardation by a 
frictional force proportional to the normal reaction is precluded. The 
analysis here presented assumes that the effect of the rolling-surface in 
retarding the forward motion can be represented by a couple K, of constant 
magnitude K, acting on the bow], the axis of K being parallel to the ground 
and normal to the direction of motion. Comparison with observation has 
shown that the resulting mathematical model is adequate for a wide range 
of initial and rolling-surface conditions. In particular there is close agree- 
ment between theoretical and observed behaviour in the game of bowls. 
Some details of the comparison are given in section 10. 


2. Frames of reference 

Fig. 1 shows the axial frame Oxyz, fixed on the horizontal rolling-surface. 
O is the initial point of contact of the bow] with the surface, Oz is the initial 
direction of motion and Oz is vertically upwards. The distance OP 
travelled along the path in time ¢ is s, the range of the bow] is R, and the 
total angle of deflexion A. The coordinates of the centre of mass @ of the 
bowl at time ¢ are (x,y,z). 

A vertical section of the bowl is shown in Fig. 2. The moving axes Génf, 


[Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 3, 1958] 
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have Gé in the direction of motion, Gy along the axis of symmetry, GZ in 
the vertical plane containing Gy. A, B, A are the corresponding principal 
moments of inertia. The angular orientation of the bow] at time t is specified 
by the Eulerian angles 0, %, 6 where @ is the inclination of Gf to the vertical, 
ys is the angle between the plane (nf and its initial position, ¢ is the angle 
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Fic. 1. Bowl trajectory and fixed frame Oxryz 


between the plane Gy and a plane fixed in the bowl containing Gy. These 
angles are depicted on the unit sphere in the auxiliary diagram in Fig. 2 
By analogy with a spinning top (Deimel (2)), the angular velocity d:b/dt is 
called the precession. From measurements of the curvature, the curve DPE 





Fic. 2. Bowl section and moving frame Génf, and Eulerian angular 
coordinates 0, pb, 


is very closely the arc of a circle. Let the centre be Q; let GQ = c, and 
LDQE = @. 

Fig. 2 shows a point contact P between the bowl and the rolling-surface, 
giving an overturning moment equal to Mgcsin(@—@).. Since the bow! 
forms 4 small depression in a deformable rolling-surface this moment wil! 
be effectively reduced throughout the motion. The reduction may be 
adequately represented by ascribing toc a value less than its actual value. 


In the present paper this is the only way in which effects of a deformable 
rolling-surface will be considered. 
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The six coordinates x, y, z, 0, %, ¢ are connected by three constraint 
relations arising from the conditions of no slipping at P. It is easily verified 
that one of these relations is integrable, enabling z to be expressed in terms 
of @. The other two constraint relations are non-integrable, so that the 
dynamical system is non-holonomic and requires the remaining five co- 
ordinates to specify its configuration. 

For a typical 5-inch bowl, measured data are 
M = 3-30 lb,a = 0-208 ft, A = 0-0542 lb ft?, B = 0-0586 lb ft?, @ = 15° 


> 


where is the mass and a the maximum radius. 


3. Dynamical conditions 
In the following analysis dots indicate time derivatives. If u is the 
velocity of the bowl when P is at O, we take as initial conditions 
t= 0, ex z=y= 0, == 4= 0, 
s=ag=t=u, y=0, O=Y=0. 
At the end of the trajectory all variables are denoted by the corresponding 
capital letters. The end conditions are 


t= T. = 8, ~= X, y= Re 6= 0, m = ¥, 
Observations of the motion of bowls provide important information on 


the relative magnitudes of the angular coordinates and their derivatives. 
This information may be summarized as follows: 


# remains very small until close to the end of the motion, 

‘Y may exceed 47, 

6/% is small throughout the motion, 

0, i <¢ except close to the end. 
Thus, at appropriate stages of the analysis, terms involving products and 
powers of @ and & will be neglected and approximations with respect to 0 
made. A solution which shows the general character of the m>tion can, 
of course, be achieved by treating @ as a first order small quantity at the 
outset. This solution, however, has serious errors in the end conditions. 


4. Kinematical relations 
In the following, the components and rates of change denoted by a dot 
are taken relative to the moving frame Géyf. Thus from Fig. 2, the angular 
velocity w of the frame Géf has components given by 
mee idk 
—wsin 0 
weos 6 


092 .43 Au 


(4.1) 
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Also, the angular velocity Q of the bowl is 


Q = $1. 
d—wsin 0 (4.2) 
w& cos 6 
Let c, = ccos 0, c, = csin 0, so that from Fig. 2, writing r = GP, we have 
rx 0 
(a+c¢,) sin @—c, (4.3) 


—(a+c,) cos 0+-c¢, 


For brevity we introduce 


f,i(@) = (a+¢,) cos 0—e,, (4.44) 
f2(9) = (a+¢,) sin 6—cy, (4.4b) 
f,(@) = a+e,—c cos(O—@) 

= f,(0) cos 0+-f,() sin 8, (4.4¢) 
S,(@) = a+c¢,—2c cos(O—86) 

f,(¢)—c¢ cos(O—8@), (4.4d) 

Thus, we have 

f,(@) = —(a+e,)é sin 8, (4.5a) 
f,(8) = (a+-c,)0 cos 8, (4.5b) 
f,(0) = —césin(Q—8). (4.5¢) 


By (4.4a) and (4.4b) we can express (4.3) as 
r= 0 . 
f.() (4.6) 
—fi() 
For use in the no slipping conditions we shall require the components, 


in Gn, of the velocity of P relative to G, ie. QAr. From (4.2) and 
(4.6), on simplification we have 


QAr = [—f,(0)p+epsin(O—8)]. 
—f(9)0 
—f2(0)0 
It follows that 
(@Ar) i ~f,(0)b+(a+-c,)66 sin 6 +e sin(O—0)—cy6 cos(@— 8) 


—f1(0)8-+-(a+-¢,)6? sin 
—f(0)8—(a-4 c,)0? cos 0 


Also, from (4.1) and (4.6) we have 


o.F= —fx(O)w. (4.8) 








THE DYNAMICS OF A BOWL 355 


5. Equations of motion 

One of us (M. N. B.) has derived the equations of motion using the form 
of Lagrange’s equations appropriate to the non-holonomic case. For 
brevity here we obtain the angular equations of m tion directly by the 
usual moving axes method. 

Let the operation d'/dt denote a rate of change relative to the fixed 
frame G£'n'C’ which is instantaneously coincident with the moving frame 
Génl. Let 

h = angular momentum of the bow] about G, 

v = velocity of centre of mass G, 

F = force on the bow] at P, 

K = retarding couple on the bow] as in section 1. 


It is easily seen from (4.2) and Fig. 2 that 


h = —Aé ‘ K => 0 . 
B(¢—wsin 8) —K cos 6 
Ayscos 6 —K sin@ 
and that the acceleration of gravity ¢ is 
g = ae (5.1) 
gsin 6 
—g cos 6 


The condition of no slipping at P is 
v+Qar= 0. (5.2) 


The linear equation of motion of the bowl, of mass 1, is 


mY — Fig, (5.3) 
dt 
d'y dv 
‘he —— = —-+ wv. 
where di di oe 
Also the angular equation of motion of the bowl about its centre of mass 
is : 
d h —-rAF-+K, (5.4) 
dt 
é lh 
where = —_ a +wAh. 


Eliminating v from (5.3) by (5.2), we have 


u| i “(2 r)—w A(2 \0)| — F+ Mg. 
d 
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Hence it follows that 


rAF = —Mr/ Re Ar)+wA(Q n)+8| 


= —Mr/ Re \r)+(w.r)Q-4 e|. (5.5) 
Substituting (5.5) in (5.4), we have 
amare : h+Mr A Fe \r)+(w.r)Q 4 ¢|—K == §), (5.6) 
€ C 


Now, using the component forms (4.1), (4.2), (4.6), (4.7), (4.8), (5.1), and 


performing the remaining differentiation and vector operations, the three 
component equations are 


{Bcos 0+ M f,(9) f4(8)}ox — Mge sin(@—6)— 
—{(B—A)sin 6 cos 0+ Mef,(0) sin(Q@—6)\y2+ 
+{4+- Mc?+ M(a+e,) f,(0)}0—Mc(a+c,)@ sin(O—@) = 0, (5.7) 
(B+ M[f,(8)"}b—M(a+c,) f,(0)60 sin 0@— 
—{Bsin 0+ Mcf,(4) sin(O—@)b— 

—{Bcos 6+ Mf,(@) f,(8)}0+ K cos@ = 0, (5.8) 

Mf,(0) fo(0)6—{ B+ M(a+-c,) f.(8) sin 0}66+ 

+-{A cos 0— Mef,,(8) sin(O—6)s— 
—{(24—B)sin 0+ Mf,(0) f,(0)\W~O+K sind = 0. (5.9) 


For convenience we obtain the angular equation of motion about the 


vertical through P by subtraction of sin@ times (5.8) from cos@ times 
(5.9). This yields 


—{Bsin 0+ Mef,(8) sin(O—8)}d— 
—{Bcos @—Mc(a+c,) sin(@—8) sin 0}66-+ 
+{A+(B—A)sin?é+ Mc? sin?(@—6) b+ 
+-{2(B—A) sin 6 cos 0+ Mcf,(8) sin(O—6)}yd = 0. (5.10) 


In the next sections approximate first integrals which are sufficiently 


accurate to give a satisfactory solution of the problem are obtained from 
(5.7), (5.8), and (5.10). 


6. The principle of angular momentum about a vertical axis 
For a typical bowl both c and B—A are very small. Since sin @ and 
sin(@—@) are also small, an approximate integral of (5.10) is 
— Bésin @+-Ay’ = constant, 


i.e. the angular momentum of the bowl about the vertical through P is : 





veins tt boaabi oS ca A boon aR 
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constant under the adopted assumptions. (The ‘rolling-surface moment’ 
about this axis is of unknown form, but is certainly very small, and has 
been neglected as indicated in section 2.) Using the initial conditions 
0 = 0, » = 0 the equation becomes 


Ay = Bosiné. (6.1) 


7. A momentum-impulse equation 
The derivation of an approximate integral of (5.8) is facilitated by 
rewriting it in a different form using (6.1). Since A/B = 1, we take here 
b = osiné 
so that M(a+e,)f,(0)b0 = M(a+c,) f,(0)66 sin 0. 
Using this, (5.8) may be written as 
| B+ M[f,(8)?}$—2M (a+e,)f,(0)G6 sin 6— 
—{Bsin 0+ Mcf,(6) sin(Q@—6)b— 
—{Bcos0—2Mef,(0) cos(O—@)}y~b+ K cos? = 0. (7.1) 
On neglecting terms in ¢ and integrating we have 
t 


{B+ M{ f,(0)}?}6— By sin 0+K [ cos @ dt = constant. (7.2) 


0 


It is satisfactory to simplify this further by writing cos@ = 1. Thus 


/,\(@) = a, and using the initial conditions & = uja, y = 0, (7.2) becomes 


(B+ Ma?)¢— Bysin 0+ Kt = (B+ Ma*)u/a. (7.3) 


This is an approximate form of the momentum-impulse equation for the 
axis through P parallel to Gy. 


8. The principle of energy 

An energy integral may be obtained from (5.7), (5.8), (5.10) in the usual 
way. We multiply these equations by 6, ¢, w respectively, add and obtain 
on integration 


{B+ M[ f,(0)}*}62—{B sin 6+ Mef,(6) sin(O—8)}op+ 


t 


+K | 6cos0 dt +4{A+(B—A) sin’ + Mc? sin*(O—6)\y?+ 
0 


+}{A+Me?+ M(a+e,) f,(0)}6?—Mge cos(O—#) = constant. (8.1) 


The equation (8.1) is the principle of energy, and includes the work done 
in time ¢ against the retarding couple K. An excellent approximation 1s 
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obtained by neglecting the very small terms in y? and 62. Since by (6.1) the 
term in ¢y is also eT to the same order, (8.1) becomes 


1{ B+ M{f,(0)?\¢?+-K l cos 6 dt —Mgccos(Q—8@) = constant. (8.2) 
0 
This approximation neglects the kinetic energy of tilting and precession. 
As in section 7 we put cos @ = | in the first two terms of (8.2), and using 
the initial conditions 6 = @ = 0, 6 = u/a, we have 


1(B+ Ma*)¢?+ K¢—Mge cos(O—0) = 4(B+ Ma®)u?/a?—Mge,. (8.3) 


9. Solution of the differential equations 
On writing B’ = B+ Ma? the first integrals (6.1), (7.3), (8.3) are expres- 


sible as B. 
y= qesin 6, (6.1 a) 
- ww K B . 
= ———t4 — gain, 7.3: 
, a B an satis 
™ u= eK, Mg : , Boal 
1g? = a Bet pr (008(9—8)—e]. (8.3) 


On transposing (B/B’)wsin@ in (7.3a), squaring and neglecting terms 
in y* (by virtue of (6.1a)) we have 
2 2 
jen Se Fp. (9.1) 
a® aB B? 
while integrating (7.3) over the interval (0, t) gives 
u K 


t 
B 
d= ab Ta P+ = a [isin e at. (9.2) 


thi 


9.2) in (8.3a) it reduces to 


_ 


On substitution for ¢? and ¢ from (9.1), 
t 


_ MgB’ 
[ wsin @ de = FR [c cos(Q—6) —c,]. 


0 


Differentiation with respect to ¢ gives 


dys 1 sin(@—8) (9.3) 
dd p @sind ’ this 
KB 
where EE seo 9 
P= BMg0c sate 


Integrating (9.3), after writing sin(Q—0) = @—8, sind = 0, we have 


l 
= = (log 6— 3) + constant. (9.5) 
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Since 6/0 < \log@|, except close to the end of the motion, and 6/pOy is 
found to be small, it is satisfactory to write (9.5) as 


1 
% = —log@é+-constant, 
ha 


whence, using the end conditions 4 = VY’, 6 = ©, we have 
0 = Oerd-*), (9.6) 
This equation does not satisfy exactly the initial conditions 6 = 4 = 0 but 
#/© = e-”¥ would be expected to be very smali. This is verified by observa- 
tions of p and ¥’ given in section 10. 
We next consider the first integral (6.l1a). Using the approximation 
c = 0, cos@ = 1 we have 


$ = ad. (9.7) 
Using (9.7) and putting sin@ = @, (6.1a) reduces to 
o et (9.8) 
di Bae 


Thus, approximately, the radius of curvature of the trajectory varies inversely 
as the angle of tilt of the bowl. 
Equation (9.8) becomes, on substitution for @ from (9.6), 


ds _ aA ne_y, (9.9) 
dp BO 
The solution of (9.9), using the initial conditions s = 0, / = 0, is 
s= Lv —e-Py), (9.10) 
q 
where _ PBO_ _KB (9.11) 


~ ad aAB'Mge 


Equation (9.10) is the intrinsic equation of the bowl path. Expressions for 
the coordinates 2 and y of the centre of mass suitable for computation of 
the bowl trajectory can now be obtained. For convenience we introduce A 
defined by A= ery, (9.12) 

so that the range of A corresponding to 0 <s < S is, from (9.10), 

From Fig. 1, ¢ = écosy%, y = sing, so that 
t wb 
atiy = [ sew dt = [ elt 
n) ny 
i s ert — 
= | wa Bi tll {1—A(cos¥+-isin })}, 
q(p?+ 1) 


ds aAert 
ay! = BO” 





w 
| civ db, by (9.9), 
0 











360 M. N. BREARLEY AND B. A. BOLT 


on substituting from (9.11), (9.12). Equating real and imaginary parts we 
have 


Pp epY p ) 
r= ron ow {p—A(pcosy—sin #)}, (9.13) 
i 
y= ai — {1—A(cos /+-psin y)}. (9.14) 


For computation it is convenient to introduce the range R = ,(X?+-Y?) 
into (9.13), (9.14). Now from (9.13), (9.14), inserting the end conditions 
a= X,y=Y,\A =e", we find 

7 p ert 


‘3 Vip+1) q’ 


omitting second-order terms. Thus (9.13), (9.14) may be written 





== <a P—Npeosy—sin f)}, (9.13a) 
V(p?-+ 
m5 — 
y= (pPpay tt —Alcosd-tp sin )}, (9.14 a) 
V(p?+ 


where the angle y& is, from (9.12), given by 


b= Flog, (9.15) 


—-1 


for l1>A> [Avot +008] 


Hence for particular values of p and R values of coordinates (x,y) corre- 
sponding to a suitable set of values of A may be calculated. In particular, 


since A is very small close to the end of the trajectory, we have with 
sufficient accuracy 


fe! TS i Foe Te | 9.16) 
v(p?+1) v(p?+1) : 

> 1 Os . 

¥ = Flog| w+ DRI. 9.17) 
p lp? ( 


10. Discussion 


The formulae giving x, y, and % for a particular bowl depend upon the 
range # and, through the parameters p and q, the moment K. Since we 
would expect K to vary with surface conditions it must be estimated for 
any given surface as follows. From (7.3a) and (9.7), neglecting the small 
term containing ¥sin @ we have 


é= u——t. (10.1) 
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Since A has been taken as constant the acceleration of G is thus approxi- 
mately —Ka/B’. Thus 

Ka 2B’ 8 

——aee? ena K =... 
2B’ a T 
In practice it is sufficiently accurate to replace S by R in (10.2) since 
experiments on various types of surface give 0-9 < R/S <1. Table 1 
shows values of R/T? obtained with the same bowl on three bowling green 
surfaces described as ‘fast’, ‘medium’, and ‘slow’ respectively. 


| (10.2) 


TABLE | 
Observed ranges (ft) and times (sec) 











Trial Green 1 (fast) | Green 2 (medium) Green 3 (slow) 

number) RK | T | R/T?| R rT |r| R | 7 [RP 

g°2 45 | O45 | 11's 4°9 | 0°48 | 13°71 4°6 | 0-62 

2 174 | 6&3 | O44] 163 | 5:7 | O50] 219 | 5:7 | 067 

3 3°°7 86 | 0-42 | 263 72 | O51 | 29°6 67 | 0°66 

4 47°7 | 102 | 0°46 | 3571 8-2 | 0-52 | 37°6 76 | 065 

5 57°7 | 11°7 | O41 | 43°0 88 | o56 | 46°6 gt | 056 

0 64°74 | 12°5 | Og | 54°3 | 103 | O51 | 554 99 | 057 

7 78-2 | 13°6 | o-42 | 61-9 | 10°9 | O52 | 65:0 | 105 | O59 

5 89°4 | 14°7 | o-41 | 72-2 | 11-9 | O51 | 74°8 | I1-O | O62 

9 97°3 | 14°9 | 0-44 | 83-4 | 12°3 | O55 | 81-0 | II | 0°66 

10 gt-2 | 12°7 | 0-56 | gor5 | 12:1 | 0°62 
II 100°6 | 13°6 | 0°54 
































The values of R/T?, by (10.2), indicate that K for a particular green and 
bowl is essentially independent of R and hence of the initial speed of the 
bowl, at least for the ranges investigated. This provides experimental 
support for the hypothesis introduced in section 1 that K may be taken 
constant for the same bowl and bowling surface. From (10.2) and the mean 
values of R/T? given in Table 1, the values of K for each green are 0-84, 
Ol, and 1-20 ft pdls respectively. 
Insertion of the end conditions = 0, t = 7 in (10.1) shows the total 
time of travel of the bowl to be 
Bu 
= 
Ka 
On substituting this in (10.2) and writing S as R, the range achieved for a 
speed u of projection is found to be 
Bt u. 
2Ka 
Por given u, the range is therefore directly proportional to the speed of the 
green as measured by the factor 1/K. 
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To plot trajectories for a particular bowl and bowling surface an appro- 
priate value of p must be found. This may be determined as follows. The 
angle of deflexion A, shown in Fig. 1, is given by tanA = Y/X, in which 
substitution from (9.16) gives 


tanA = p=. (10.3) 


Thus p may be found by measurement of A in a trial trajectory. 

If the value of p so found is used in (9.4), this equation then gives the 
value of c appropriate to the bowl and surface concerned. We note that 
for a 5-inch bowl numerous dynamical experiments on various bowling 
greens show that if the values © = 15°, c = 0-0016 ft are used in (9.4) to 
obtain p as a function of K, satisfactory agreement between calculated 
and observed paths is obtained. In agreement with section 2, this value 
of c for a deformable surface is considerably smaller than that obtained 
from statical experiments on a polished steel surface, namely 


c = 0-01+0-002 ft. 


The equation (10.3) shows that for the same bowl and rolling-surface th: 
angle of deflexion is the same for all ranges R. 

This theoretical prediction was tested by measuring values of (X, Y) for 
each range on the three greens used in Table 1. When plotted the values 
showed no systematic trend from a linear relation. The least-square 
solutions with their standard errors were: 


Green 1. Y = (0-182+0-003)X —(0-57+0-21) 
Green 2. Y = (0-156+0-004)X + (0-60+-0-20) 
Green 3. Y = (0-141+0-003)X —(1-51+.0-33) 


The result provides a very satisfactory check on (10.3) and hence on the 
approximations of sections 6-8. The corresponding values of A are 
10-3°.0-1, 8-9°+0-2, 8-0°+0-1. The decrease in the angle of deflexion 
with the decrease of speed of the green, by (10.3), entails an increase in p. 
This is in agreement with (9.4). 

Finally, we give in Fig. 3 a comparison between the calculated and 
observed trajectories for two indoor test-tables. 

The test-tables, one at Sydney and the other at Adelaide, are canvas, 
over sponge rubber, over slate. They are used to check the bias of bowls 
used in competition against the bias of a standard bowl. Fig. 3 shows 
excellent agreement between theoretical and observed trajectories in both 
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cases. In particular, measurements of ¥’ at the end of each trajectory are 
of the right order, as shown in the following table. 






































R A K p c YW calc. | ¥ obs. 
(fi) | (degrees)| (ft pdls) (ft) |(degrees)\(degrees) 
Sydney table 252 22°5 0°397 2°41 | 00018 106 100 
Adelaide table 315 22°6 0°392 2°40 | O-0017 113 119 
¥4 
SS 39s asec Sydney test table 
vr and master bowl. \ coicuiated paths 
Adelaide test table 
10+ and master bow!. j « 
+ ™% Observed co-ordinates 4 
sr ta 
b} 
4; 
2b 
' 
nk — 








l —— | i 
4 é oo 2... a ee 


Fic. 3. Comparison of calculated and observed paths 


These results show that the treatment presenved is adequate to deal with 
the dynamics of a bowl under normal conditions. 
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SUMMARY 

Methods are developed for representing numbers in a modified ternary scale, 
which enables the operation of multiplication to be fabricated in an automatic 
binary computer with the minimum number of additions or subtractions. Two cases 
arise according to whether a parallel or serial machine is being used. In the former 
case, the number of additions (or subtractions) can be reduced to an average of 
a third the number of digits in the multiplier with a maximum of a half the number 
of digits. In the serial case a multiplication by a 2p-digit multiplier can be achieved 
as a single operation with only p adders. 

The application of this technique to division is also considered. 


1. Introduction 

In both serial and parallel automatic computing machines, the operations 
taking the longest time are usually multiplication and division and these 
occur sufficiently often in programmes to make the reduction of multiplica- 
tion and division time a worthwhile project. 

In fact many schools of programming at present estimate the time of 
any given problem by counting the number of multiplications and multi- 
plying by some fixed multiple of the multiplication time. 

Many machines now available or under construction already have 
additional mechanisms to speed multiplication. 

This paper gives a survey of some possible schemes for fast multiplica- 
tion and gives the necessary mathematical analysis to effect further im- 
provements over those already incorporated or envisaged for the present 
generation of automatic computers. 

This analysis hinges on the representation of a binary number as the 
sums and differences of a set of binary powers. 

From this analysis it is possible to derive optimum schemes of division. 


2. Multipliers in serial machines 
The simple multipliers in serial machines are based on the binary 
analogue of the continued addition method used in desk calculators. 
This consists, for each digit in turn, of either (i) adding the multiplicand 
to the partial product and then shifting the latter or (ii) just shifting the 


partial product, obeying (i) or (ii) according to the value of a digit of the 
multiplier. 


[Quart. Journ. Mech. and Applied Math., Vol. XI, Pt. 3, 1958] 
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The inspection of each digit of the multiplier is usually achieved by 
holding this in a circulating register n+1 digits in length (n-digit 
numbers). In each cycle the binary pattern is shifted back one place and 
by inspection of a fixed digital position at a fixed time in each cycle the 
digits are inspected in turn, the least significant digit first. 

Similarly, the partial product can be shifted back one place per cycle. 
An adder is provided in its circulation path to allow the multiplicand to 
be added. The least significant digit in each cycle is deliberately lost unless 
a double precision product is required when these digits can be stored 
in the used positions of the multiplicand register. 

Thus about n cycle times are required for a multiplication independent 
of the digital pattern of the multiplier. This is a twofold improvement 
over the original design of multiplier which used registers of length 2n+-1 
to hold the partial product and the multiplier and each cycle of the 
multiplication took two machine cycles. 

To reduce the multiplication time below n cycle times two alternative 
arrangements are possible: 

(a) arrange that the shift without addition which is required if the 
multiplier digit is zero can be effected without a complete circula- 
tion cycle; 

(b) add the increments to the partial product arising from several digits 
of the multiplier in each machine cycle. 

The alternative (a) involves the provision of adjustable delays which 
can be inserted in the two circulation loops to shift these numbers over 
a number of places determined by the length of a chain of zeros in the 
multiplier. Mechanisms must be provided to sense the length of the next 
most significant chain of zeros before each addition cycle so that the delays 
can be set to the required value during that cycle. 

The alternative (b) involves a multi-input adder in the circulation loop 
of the product register, each input feed having a gated and suitably 
delayed copy of the multiplicand. The contents of the multiplier must be 
used to control the gates on these inputs. 

In practice the multi-input of say n inputs is constructed from a chain 
of n ordinary adders with the output of each feeding as one input to the 
next. Thus to deal with p digits of the multiplier in each cycle p adders 
are required. 

A second method of applying alternative (6) is to use a set of adders 
to produce a unit giving on 2” different outputs all possible 2” multiples 
of p digits. For each successive pattern of p digits in the multiplier, the 
appropriate output is gated to the input of the ordinary adder in the 
circulation path. This is an uneconomic method for large p since the unit 
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will require 2”-!—1 ordinary adders (one to produce each of the odd 
multiples) giving a total of 2?-! adders compared with p adders in the 
first scheme. 

In alternative (a) the multiplication time is directly proportional to the 
number of ones in the multiples and so assuming that all n-digit numbers 
are equally likely, the average number of cycles is $n and the multiplica- 
tion is, on average, halved in comparison to the original simple multiplier. 

Any further decreases in multiplication time by this method must 
depend on finding a representation of the multiplier which contains fewer 
ones than the normal binary representation. Since a binary adder-sub- 
tractor switchable from one function to the other at will is a well known 
device which is not much more complicated than a simple adder, it is 
natural to consider representations of the form 


n 

4 8 » > 
> (—)*d,2" (s8,,a 
r=0 


, 7, 9) 

for the multiplier which will correspond in the multiplication process to 
adding or subtracting the multiplicand (according to the value of s,) when 
d, = 1. 

For this form of multiplication unit the problem arises of choosing the 
form so that the number of non zero d’s is minimized. 

The multi-input adder type of serial multiplier can be generalized by 
replacing the adders by adder-subtractors and using a similar form for 
the multiplier to control both the gates and the function of the separate 
adder-subtractors. 

If p adder-subtractors are provided the problem arises of determining 
the range of multipliers that can be encoded with not more than p non- 
zero digits. 

To produce a workable scheme of this kind each delayed version of the 
multiplicand must be fed to a fixed input of the adder and since, if the 
scheme is to have any value, the range of multipliers will have more than 
p digits, some or all of the inputs will be fed with more than one version 
of the multiplicand. We require that these multiple inputs to each adder 
input can be mixed in a simple gate and so the problem arises of allocating 
the digital positions of the multiplier among the p inputs in sets so that 
every multiplier in the possible range can be represented with at most 
one non-zero digit in each set. 


3. Multipliers in parallel machines 

The original von Neumann parallel arithmetic unit works on the familiar 
principle of the Brunsviga machine with the exception that multiplication 
proceeds the least significant digit first and the least significant digits of the 
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product, if required, are stored in the multiplier register (1). This arrange- 
ment ensures that only a single length adder is required. 

In some machines the process of adding and shifting takes place in a 
single phase. If the addition is not required the phase will be considerably 
shorter in time since no carry propagation time is needed. 

Other machines separate the two processes of adding and shifting and 
once again the cycle time for zero multiplier digits is shorter than for 
non-zero ones. 

In either case if the shift time can be neglected compared with the 
addition time the multiplication time is roughly proportional to the 
number of non-zero digits in the multiplier. 

Subtraction is often performed by complementing the contents of the 
multiplicand register and if this also takes a time negligible to the addi- 
tion time, the possibility arises of using the form discussed in the preceding 
section for the multiplier which will further reduce the multiplication 
time. 

The possibility of achieving further decreases in multiplication time by 
using multiple input parallel adders need not be considered, since the 
cost of such units rules them out of serious consideration. 

The remainder of this paper is concerned with the mathematical analysis 
of the representations of binary numbers in the various forms and 
an investigation of the recursive type of representations suitable for 
mechanization in multipliers. 

The analysis is then applied to the mechanization of the division 


process, 


4. The representation of binary numbers 

It is clear that the discussion can be restricted to integer multipliers 
since the products formed using other multipliers can be derived from 
that of a corresponding integer by a simple shift. 

Suppose a given integer N is represented as 


or > 62 (b, = 0,1; r = 0,1,...). (1) 


For n-digit numbers 6, = 0 (r > n), but in the following analysis which 
gives recursive relations the finite nature of N is irrelevant. 
We require to find the condition on s, and d, such that 


vt - a Aaa > (-" d, 2" (b,,8,,d, = 0,1; 7 =0,1,...). (2) 
From (2) we have 
> {(—)* d,—6,}2” = bo—(—)* dy = (say) (3) 


r=1 
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where e, can take the values —1, 0,1, 2. But 
€é) = 9 (mod 2), 
i.e. b, =d, (mod 2), 


and since b,, dy are binary elements, 


Now define gfe) = 0 (e #2, i.e. e = 0,+1) | 5) 
and g(a) = 1 (otherwise) I" 
Then (3) becomes > {(—)* d,—b,}2"1 = g(eq). (6) 
r=1 
Hence > {(—)*d,—b,}2"-1 = g(eo)+b,—(—)"d, =e, say. (7) 
As before, e; =90 (mod 2) 
d, = {b,+ (e)} (mod 2). (8) 


This shows that e, can only take values 0, 2, and applying the above 
process successively we obtain 


S {(—)d,—b,}2" = f,+b,—(—)"d, =e, say, (9) 
where re fi = 9-1), 
whence d, = {b,4f,} (mod 2). (10) 
Since e, = 0, 2 are the only possible cases it is easy to verify that 
Seca = fb +84). (11) 
Equation (10) may be written in the forms 
f, = (b+d,) (mod 2) = b,+d,—2b,d, = (b,—d,)?, (12) 
d, = f,+-b6,—2f,b, = (b,—f,)? (13) 


since x? = x for a binary element. We can eliminate f, from these equa- 
tions. Substituting (12) in (11) 


Sisx = (by +-d,— 2b, db +-8,4, 
= b,—b,d,+-8,d,. 
Replacing t by ¢+-1 in (13) and substituting for f,,,, 
d,., = (b,—b,d,4+-8,d,)(1—2b,,,)+5,,, 
= by +b, ,., —2b,b,., —(b,—8,)(1 — 2b,.,)d, 
= (bb 444)? — (bj —8,)(1 — 2b, 1). (14) 





Thus for any arbitrary sequence of s,, we may represent any integer .V 
given by (1) in the form (2) by a sequence d, satisfying d, = b, and the 
recurrence relation (14). 
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5. The minimal digit representation 


We require to choose the sequence s, so that ¥ d, is minimized for any 
integer V, where 


Sd, =dot ¥ dpy = bot ¥ (6,—b,.,)— ¥F (b,—8,)(1—25,.4)d,. 
r=0 r=0 r=0 r=0 

For a fixed n-digit N, b, = 0 (r > n) and the first two terms are a finite 

fixed sum. Hence ¥ d, is minimized if > (b,—s,)(1—2b,,,)d, is maximized. 
This in turn is maximized if c, = (b,—s,)(1—2b,,,) is maximized for 

each non-zero d,. c, = 1 if and only if 


(i) KX =1, = 0, b 41 = 0 


r r 


or Gi) 6,=0, 4=1, 6,,=1, 
i.e. 6, + b.,, = 4, 
if 6, = 6..,, 


c, = (b,,,—8,)(1—2b,.,) = —(b,.,—8,)?. 
In this case ¢, is maximized at value 0 if s, = b,,,. 
Thus the sequence s, = 6,,, gives the minimum value of } d, since 
this maximizes c, for all r and hence for all r with d, + 0. 
For these values of s,, (14) leads to 


dog, = (6,4, —5,)?+ (6,4 —6,)(1 — 2, .)d, 
= (b,,,—,)?(1—d,) (15) 
8, = _ (16) 
Thus, the minimal digit representation D is given by the rule that a non- 
zero digit occurs in D if and only if the corresponding digit in the binary 
representation differs from its lower neighbour and the preceding digit 
of D is zero. The digit is positive or negative according to the value 0 or 1 
of the higher neighbouring digit in the binary representation. 
A direct proof that this rule gives a minimal digit representation can 
also be given.t ; 
Denote the sequence of d’s given by the rule s, = 6,,, by {d,} and the 
sequence by any other rule {d,}. 


Let T =~2¢ 


n 
Then Ty = 0. 
It is required to prove that 7), > 0 for all n such that ,.1 = 6, and 
for all sequences {s,}, {b,}. 


+ This alternative proof was added to the original paper at the suggestion of a referee, 
who was not convinced by the earlier argument. 
5092.43 Bb 
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Consider the following transition table: 








| ‘ min(d},,, @}') ,)— 
d, d, uma! (8, = bpsa)| Ops (8 F bps) dy. dys 
0 . | Maar (Bpsy—by)* (bpsy—by)® 0 
0 l | (by. <b, )* (bps, —bp)? 0 (bps1- by)? 
1 0 | 1 (bp, —5,)? — (bp, —0,)? 
1 1 ; 1 0 0 














Thus 7,,,, > 7, unless d, = 1, d,,, = Oand b,,, # 6,. Suppose nj, ng... 
are the ascending sequence of n’s for which d, = 1,d,,,, = Oandb,,,, + b,. 


Then T+ = T,,—1. 

Now T,, = T,,-1+1, 

thus T,,+1 = F,-1° 

Hence T,2>T%=0 (n = 0,1,..., my—1), 
T, > Tr = T1290 (wn = 2, 4+1,...,m.—1), 


i ST. =7,.,290 (n= 2,+1,...,%4,;—1). 


If n is not in the set n,, ng,..., 


If b,,., = 6,, n is not in the set n,, g,.. 
Therefore, if 6, ., = 6,, T,, > 0. 

This rule is particularly easy to mechanize since it only involves scanning 
three digits of the multiplier and storing the preceding d digit. 

Since no two successive d digits can both be non-zero, in a parallel unit 
each actual addition can be followed by a double shift and the d-digit 
store eliminated. Similarly, in the analogous serial unit the variable shift 
can be adjusted to give the same result. 


n> 


6. The number of digits in minimal digit representations 
For any given N, > d, can be calculated from 


> 4, = b9+ Dd 4, 44 = bo+ > (6,41—6,)%(1—d,) 
= byt ES Opa—b,)?— ¥ (bes —b,)*4, 
= bot ¥ (bss —b,)?— ¥ (bps. —5,)%(6, api 1)?+ 
+3 (b 2(b,—b,_,)*d,-, 
= bot > ( P— >  )* 2450) °( : je 
~ )*( )®( )?( P+ 
Now > (6,,,—6,)? is the number of changes of value in the sequence 


bo, by,..+ > (b,.,—5,)*(6,—b,_,)® is the number of chains of alternating terms 
of length 3, > (6,.,—5,)?...(b,.3—0,42)? the number of chains of length r—s. 
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Hence the number of non-zero digits in the minimal digit representa- 
tion of V is the number of chains of alternating terms of even length less 
the number of chains of alternating terms of odd length. 

For example if VN = 1466 = 10110111010 we have 





uw 


Chain length . , 2 3 4 
Number of chains. s 5 2 0 

















The number of digits required is 8—5+2—0 = 5. The coding proceeds 


as follows: 





’ Oo I 2 3 4 5 6 | 8 9 10 II 12 13 

¢ I ° I I I ° I I ° I ° ° ° 
Ab)? I I I ° ° I I ° I I I ° ° 
i I fe) I ° ° I ° ° I ° I ° ° 

I ° I I I ° I I ° I ° ° ° ° 


N = 2u_ 99269342 
= 2048—512—64—8+42 = 1466. 


The chain rule is rather cumbersome to operate and a count in the 
recoded form is in practice a quicker method. 

For any assessment of the efficiency of the method, however, it is the 
average number of digits occurring which is required. If it is assumed 
that all m-digit numbers are equally likely the result may be obtained frem 
the total number of »..n-zero digits in the representations of all the in- 
tegers in the range [0, 2"—1]. 

Let S,(b,d,d’) denote the set of sequences of b's for which 


(i) b, 0 (r > n), (ii) b,, — b, (iii) d,, — d, (iv) do. d' 
and let S,(X) denote the number of sequences in S,(X), where X denotes 
iny ordered triple (b,d,d’). 
From the relation d, ., = (6,.4—6,)?(1—d,), 
S(0,0,1) = S,,(0, 1,1) = S,,(1, 0, 0) = §,(1, 1,1) - 0. 


Consider the relation between S,(X) and S,, ,,(X). We have the following 


transition table: 

















bast ° Omer . 
] 7 
Case b, dy dast 1, 1 dyse Inet “n-2 
i 

A ° ° ° ° ° I . 

B ° I ° ° ° ° . 

C I ° I I ° ° . 

D I I ° ° ° | ° ‘ 
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Thus we have 


Sn4i(A) = S,(A)+S,(B)+S,(D) 

Sn+i(B) = S,(C) (17) 
Sra(C) = S,(B)+8,(C)+8,(D) [- 

S,,4,;(D) = 8,(A) 


This is a system of linear difference equations whose characteristic equa- 
tion is easily calculated as 


E(E—1)(E+1)(E—2) = 0. 
Thus S,(X) = A,2"+B,+C,(—)", where A,, B,, 
the sequel it transpires that we only require S,(€), 


S,(C) = So(1, 0, 1) = 0 


C, are constants. In 


S,(C) — i, 
S,(C) — z 


Substituting these boundary conditions gives 


S,,(1,0,1) = 32"—3—}(—)". (18) 
If d,(X) denotes the total number of non-zero d’s in the minimal digit 
representation of the sequences of S,(X), reference to the table shows that 


d,43(A) _— d,,(A)+d,(B) +d,,(D) 
d,,.,(B) _ d,(C) 


dy.x(C) = dy(B)-+d,(C)+d,(D)+8,(B)+8,(D) - 
d,,,(D) = d,(A)+S,(A) 
If d,, = 34,(2) (X) is the total number of non-zero d’s in > S,,( 
d,,.4 = 2d, +8,(A)+S,(B)+S,(D). 
Clearly S,(A)+S,(B)+8,(C)+8,(D) = 2", 
whence (EH—2)d, = 241-8, (C) = $2"414-44-}(—)*. 
The general solution of this difference equation is 
d, = A2"+ Jn2"+1_4— 1 (_)n, (20) 
The boundary condition is d, = 1 whence 
d,, = 3{(3n+7)2"+2—9+4 (—)n+1], (21) 


This gives the number of non-zero d’s in the set of n+1 digit numbers 
and must be compared with (n+1)2" for the binary set. For large n, 


d,,/(n-+-1)2" ~ 3 showing a 33 per cent. saving of digits and corresponding . 
decrease in multiplication time. 
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7. The minimum number of digits required to represent a range 
of integers 


We now consider the multi-input form of serial multiplier. Given n 
inputs we require to find N so that all integers in the range 0 < x < N 
can be represented so that $d <n. This is equivalent to determining 
the least N+-1 for which ¥ d > n for all representations. 

Consider the generating function 


G(x,y) = [J (y+14+a%y)= DS D> aay’. (22) 
i=0 ; r=—@ s=0 


Then a,, represents the number of ways of representing r with s non-zero 
elements. Write 


x 


> a,,y° = a,(y) = y™ > buy (b9 # 9), (23) 


s=0 
so that n, is the lowest power of y in a,(y), i.e. the minimal digit representa- 
tion of r has n, non-zero digits. Now 
G(x,y) = (xy+1+2x-1y)G(x?, y) ; 


@ r 


hence a,(y)a? = (xy+x4y) SY alye*+ Y a,(y)e* 


-. 


= > xaY)t+aayjer"+ SY aly. (24) 


ree @ pat 
Comparing coefficients of x” 
a,,(y) = 4,(y), (25) 
A s(¥) = Y{a,(y)+4,.1(y)}, (26) 
and comparing leading powers of these polynomials, we find 
Ne, = Ny, (27) 
Nops, = Min(n,,n,,,)+1. (28) 


(25) expresses analytically that doubling a number in the binary scale 
does not alter the number of representations with a fixed number of non- 
zero digits. 

The first appearance of the value n occurs in the sequence n,, %, m3 
at an odd position, say 2s,+1. Then n+1 cannot appear in the sequence 
until position 2r+1, where r is least solution of n, = n,,, = n (from (28)). 
One of r, r+1 is even and the least even value of r is 48,42. Consider 
the value of n, for r = 48,+1: 


Nae, +1 = MIN(Ng,,, Noe,41) +1, 


or using (27), N4s,41 = Min(n,.,n)+1. (29) 
From (28) 89,41 = m = min(n,, M,,4,) +1. 
But n, < n, ,.,, < by definition of s,. Hence 

N,, = M4, = n—1. 
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Hence the least solution of n, = n,,, = nisr = 4s, ,, and first appearance 
of n+-1 is at 8s,,,5. 


Putting 2s,,,, = t, we have the recurrence relation 


bast = 4t,1 (30) 
with boundary condition ¢t; = 1. Hence 
tr “si g4n-1+ 4. 


Thus the largest range of integers which can all be represented with not 
more than m non-zero digits is 0 < x < $(4"—1). 

Thus all (2n—1)-digit integers can be so represented and those 2n-digit 
numbers less than 101010...10 can also be so represented. 

The above proof has not used a knowledge of the form of the minimal 
digital representation. If the fact that no two successive digits can be 
both non-zero is used a simple heuristic proof follows by observing that 
a number of 2” places will contain at most n+ 1 non-zero digits (we must 
allow the possibility of d,, = 1) in its minimal digital representation. If 
d,, = 0, the maximum value that can be represented will be a number 
s, = 0, d,, = 0, dy,,, = 1 (r = 0,1,...,.n—1). This gives the result re- 
quired. 


8. The allocation of digital positions to inputs in a multi-input 
adder-type multiplier 

As explained in section 2 for a practical utilization of the result of 
section 7 each of the 2n digital positions of the possible multipliers must 
be assigned to fixed inputs of the adder. It is now necessary to determine 
the restrictions on the division of the digital positions among the inputs 
necessary to allow each number to be represented. 

Suppose, if possible, that digits d,, d,,d,(r > s > t) are associated with a 
common input. Then for all numbers in the range [0, #(4"—1)] (m = 2n—1) 
only one of d,, d,, d, can be non-zero. .Consider the representations of the 
number 2’—2'. We have 

b,=0 (u<t) 
=1 ((<ucr) 
=QQ0 (w2>r). 

For an arbitrary sequence s,, we have equation (14) which can be sum- 
marized in the following table. 








Bess b, das 
0 0 s,d, 
0 l 1—(1—s,)d, 


l 0 l—s,d 
l 1 (1—s,)d, 
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Hence for any representation of 2”—2' we have 


4,., = 8,d, (v<t—l) (a) 

= 1—4_,d,, (v = t—1) (b) 
= (1—8,)d, (t<v<r-l) (c) (31) 

= 1—(1—s,_,)d,_, (v=r—1) (d) 

= 8,d, (v 2>r) (e) 
Now if ¢ 4 0, dy = 0, hence d, = 0 (v < t) from (31a), d, = 1 from (31 b). 
If t= 0, dy = d,= 1, and hence d, = 0 by the input rule from (31d) 


(1—s,_,)d,_, = 1, ie. d,_, = 1. Hence, for (31c), d, = 1(t<u<r—1) 
and in particular d, = 1, contrary to input rule. 

Hence at most two digits can be associated with each input and since 
there are only 2n digits to associate with n inputs each input is associated 
with exactly two digits. 

Suppose digits r and t (r > t) are associated with one input and p and q 
(p > q) with another input. We can assume without loss of generality 
that r > p. Three possibilities arise: 

(i) r>p>q>t, 
(ii) r>p>t>gq, 

(iii) r>t>p>gq. 

In case (i) consider the representations of 2"—2' given by (31) d, = 1; 


d, = 0; d, = 1; t <u < r—1 and in particular d, = d, = | contrary to 
input rule. 
In case (ii) consider the representations of the 2”—2?-+-2‘—2%, 
b,=0 (u<q) do, = 8,4, (u << q—1) (a) 
=1 q<uct) = 1—8,_,d9-4 (u = q—1) (b) 
=0 (t<u<p) = (1—8,)d, qqQ<u<t—l) (ce) 
=1 (pxmu<nr) = 1—(1l—s8,_,)d_, (u=t—1) (d) 
=0 (u>r) a 6.4, (t<u<p—l) (e) (32) 
= 1—s,_,d,_, (u = p—1) (f) 
= (1—8,)d, (p< u<r—l) @) 
= 1—(1—s,_,)d,_. (u=r—1) (h) 
= 6,4, (u>r) (i) 


As before d, = 1, so that by input rule d, = 0. From (32f) d,_, = 1 
and hence by (32e) dy = 1 (tf < u < p), and in particular d, = 1 so that 
by input rule d, = 0 whence, by (32h), d,_,. = 1 and thus by (32g), 
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d,, = 1(p<u<r—1l); in particular d, = 1 which gives a contradiction. 
Thus only case (iii) can occur. Given r, p can be chosen to have any 

value less than r and not equal tot. Hence t = r—1. 

Put r = m, then the mate to m is m—l. 

Put r = m—2, then the mate to m—2 is m—3. 

Generally 2,, 2,,, are mates p = 0, 1, 2,... (n = 1). The numbers 

— 2, or 2p-4.9t_ 99 are (for any choice of r, t, p, q) all less than 2"—1, 

and so if range of multipliers to be represented is restricted to m-digit 

numbers this pairing is still necessary. 


or 


9. The choice of s, for multi-input adder type multipliers 
We introduce the parity function 


Po=1, = Press = 1—P,- (33) 
Then the input rule can be expressed as 
9,4, d,.., = 9. (34) 
From (14) we obtain 
(1—2b,,,)p,d,d,., = p,d,{(b,.,—6,)?(1—2b,,,)+(s,—6,)} = 0 
since (l1— 2x)? = 1 for binary 2. 
This is easily reduced to 
p,d,(b,.4—8,) = 0; (35) 
if p, = 0, i.e. r = 2t+-1, (35) gives no restriction. If p, = 1, ie. r = 2t, 
(35) gives by dy = 1, fig = Bains (36) 


If d., = 0, 8, can be arbitrary. Without loss of generality it can be taken 


as in (36). Hence a a (bs, 2 —b,,)*( l —dy), (37) 


i.e. equations (25) and (26) hold for even r. 


To obtain a unique formula we must add further restrictions. Suppose 
8, is a recursive function, i.e. does not depend on the value of n. Then 


do, +2 = (bars 2—bop41)* c (82p43—Doy41)(1 ngs bo, .2)dops1- (38) 


If do, = 0, 8,,, can be arbitrary. If d,,,, = 1, (bs,,,—6.,)* = 1, i.e. 
Dos, 4 be,. Two cases arise: 

(i) bo,,, = 9, by, = 1. Ifr = n—1, dy,,, = dy, must be zero and hence 
for a recursive rule d,,,, = 0 for all cases b,,,. = 0, i.e. 8,,,(dy,) = 0 for 
all cases b,,,. = 0. Hence in this case we must have s,,,, = 0. 


(ii) by,,, = 1, b,, = 0. A similar argument shows that, if the full range 
of multipliers is to be used, 


(1—8,,,,)(1—d,,) = 1, for all cases b,,,, = 0. 
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Thus 8.,,, = Oif by... = 0. If by,,, = 1, the recursive restriction gives no 


information and we therefore write 


; ; 82rt1 = Oops 2Gors1 (39) 
where g is arbitrary. 


Substituting (37) and (39) in (38) we obtain 
doy 2 = Day s2(1—Gors1 oy +1) +5 27+41(1 —2b a, ,9)[b2p-+(1—bs,)dy,]. (40) 
This recursive scheme can be simplified by the choice of g,,., = 0 or 
Jops, = 1—dy,,, which reduces (40) to 
oy = Dap. 9t+bey.4(1—2bo,,9){ba-+(1 —bo,)ds,}- (41) 
(It is not possible to choose g,,,, so that any of by, ,o, bys, bo,, dos dap 
are eliminated from (40).) Equation (41) can be expressed in a simpler 
form by introducing an auxiliary variable c,, defined by 
Cop = (1—2bg,)do,-7-bo, (42) 
Multiplying by 1—26,, and rearranging we obtain 
dy, = (1—2b»,)co,+be,, 
and hence (1—b,,)d,, = (1—bp,)co,, 
so that using (41), we find 
(1 si 2bo,+2)Cor+2 t+ Oo 12= Bop 49t-Oop41(1— 2b oy 2){B op + ( 1 —be,)do,} 
Corse = bay41{bo,+(1—bo,)e2,} (43) 
doy 49 = boys +(1— 2D oy +2)Car2 = (Dor+2—Cor +2) 
= (bo,,o+Cop42) (mod 2). 
(43) expresses the short cut method of multiplication in the scale of 4. 
The carry into the binary pair b,,,., b.,,, takes place if and only if the 
preceding pair b,,,,, b,, is the form 11, or the form 10 with a carry into 
that pair. The resulting digit d,,,, is the sum of 6,,,, and the carry 
Cor+e- 

The digit d.,,, is given by (37) and is non-zero if and only if, after 
adding the carry into the pair b,,,5, b),,. we have the form 1, 0. 

In practical multipliers the range of multiples will be restricted to 
[0,2™—1]. In this case the recursive restriction will give no guide to the 
form of s,,, in the case b,,,, = 1, 6,, = 0 since this pattern will not occur 
in the leading position. Thus we may take a more general form for 8,,;: 


Bon41 = {1—(1— bap Oar Jo i 
which ensures 8,,,, = 0 in the case b,,,, = 0, by, =. 1. 
Since s,,,, is arbitrary if b,,,, = b,, we may take 
(44) 


Bor = Dap sa Jars 
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Hence, from (37) and (38), 
dor+2 “is (bop .2—Dops1)?—Boy4(1 — 2bg,49)(1 —bg,)(1 —Gor+1)(1 —d,,) 


sanparnnmnsnesnsans neon eel 


= (bop .2—D op)? — Bop 441 — 20 9p +2)(1 Gor +1 boy s1- (45) 

One simple case arises if we choose g,,,, = 1. Equations (44) and (45) 
reduce to ; 
ds. re (bop42—Oop41)*, (46) 

Bers1 = Dops1- (47) 


Combining (36) and (37) with (46) and (47) we have the simple formulae 
dy. sp (6,.1—5,)*(1—p,d,) 


(48) 
&, = Props t+Prsr b 


. 
Another useful case is given by go,.,; = 5y,... Equation (45) then 


becomes 
‘ doy. _ (boy 42—Dop41)*(1 —B oy 44 Dap +1); (49) 
and (44) becomes Bors, = 05,418 a+e 


giving the combined equations 
dys, = (b,4:—5,){1—(P,+P,+15,)d;} | 
8, = (Pp +P p41, bps ) 
10. The representation of negative numbers 


Signed numbers are represented by adding an extra digit in the leading 
position, zero for positive numbers and non-zero for negative numbers. 
Two representations of negative numbers are possible. If —N is repre- 
sented by > 6,2", then 

(i) in 1’s the complement system is 


(50) 





b, = 1—6,, (51) 
(ii) in 2’s the complement system is 
S bp or Sb, oF = ons, (52) 


The minimal digit representation is easily extended to these systems. 
Suppose NV has minimal digit representation > (—)*d, 2”. 
For the 1’s complement system suppose the minimal digit representa- 

tion is } (—)*d'.2". Also 

(6,41—5,)? = (b,41—5,)", 

Ca = (5,41 —5,)*( 1 —d,). (53) 
Let dy = 1—by = by = dy. 
By induction it follows that d,=d,. Also 


’ , 
>> ae oy 1—6, 44 = 1—,. 


Thus 


¥ (—)eed 27 = F (—) ed, 2F = — YF (—) nd, Fr = —N. 
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Equations (15) and (16) applied to all numbers in 1’s complements will 
therefore give the correct minimal representation if the boundary condi- 


tion dy = (b—by)? (54) 
is used, where ,) is the sign digit. Alternatively, define 
b, = bh) (r< 0) (55) 


then no boundary condition is required. 
Suppose } 6,2” represents a negative number N in 2’s complements; 
then its minimal digit representation is the same as that of 5 6,2 in 1’s 


lement where ’ ‘ . 
complement wh Yb. = Yb, 2, (56) 
A similar analysis to that of section 4 shows that 6, are defined recursively 
by b= (b,~,) 

€41 = e,(1—6,) }. (57) 
&= 1 


Hence if ¥ (—)*d, 2” is the minimal digit representation of V then 
d,., _ (b,.,—6,)*(1—d,) 
d, = 1—by = by . (58) 
8,d, = b,,,d, 


We require the following easily derived relations 


e,e,—l =e, (a) 
e,6,_, = 0 (b) 
e,b,_, = @, (ec) }. (59) 


€r414, = 0 (d) 
e,b,d, = €,b, (e)} 
After some simplification we find that 
(6,44 —6,)? = (b,.1—,)?—e, 6,(1 — 26, ,). (60) 
Hence, from (58) and (59e), 
d,4 = (b,44—6,)(1—d,)—e, b,(1 — 26, .4)(1—d,) 
= (b,4,—6,)(1—d,). (61) 
From (594d), (b}.,—,.,)"°d, = 0. 
From (58), 





{b',..(1—2b,.,)+6,..}d, = {8,(1—2b, 41) +5,4:}4, 
= (b,4,—8,)*d, = 9, 
whence ed, = b,,,4,. (62) 
Equations (61) and (62) with the boundary condition dy = 6, show that 
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all numbers in 2’s complements are correctly represented in minimal digit 
form by the rules (15) and (16). 

Thus in a parallel arithmetic unit, either complementing system can be 
used without complicating the control for negative multipliers. 

In a serial machine using a multi-input adder either of the equations 
(48) or (50) can be mechanized. 

Clearly both these schemes will continue to work for negative numbers 
in the 1’s complement system with the boundary condition (54) or the 
equivalent situation in (55). 

Unfortunately, this system is not very attractive for serial machines, 
and the behaviour of (48) and (50) applied to negative numbers in 2’s 
complements must be investigated. 

For equations (48) applied to the negative numbers ¥ b,2” (in 2’s 
complements) using (60), (59e), and (33) we find 

d,.4 a (b, 15,1 —p,d,)—e, 6, l — 2b, ri( l —p,d,) 
= (b,,,—,)*(1—p,d,)—p,., e, 6,(1—26,,,). 
The second term in this expression will not vanish in the cases when 
Pry = 1, 6, = 1, 6, = 0 (8 <r). Hence the scheme (48) is not very 
suitable for negative numbers in 2’s complements. 

Similarly, it can be shown that the scheme (50) is not suitable for 
negative numbers in 2’s complements. 

We may investigate the choice of g, in equation (44) which will give a 
scheme independent of sign of number in 2’s complements. Equations 
(37) and (45) can be combined to give 

d,.4 = (6,41 —6,)?(1—p,d,)— Py +1 b,(1 —2b,.,)(1 —g,)d, | 
— Pr Opsa + Prat OG, 

For a negative number in 2’s complements (63) with (57) gives, after 

some simplification 
dy = (6, +1 —5;)*( l P,4,)—Pysr bi 1 — 2b,.1)( 1 —g,)d, 
ae (b,.4—5,)?( 1 —Pr d,)—p, +1 b,( 1 — 2b, .1)( 1 —9,)d,— 

—Pp, a(l — 2b,.,)[¢{6,+(1 ais 2b,)( 1 —g;)}+6,(g,—9;) |d, (64) 

and 8,d, — (DP, O41 t+PrirO,9;) d, 
= (Pp Opa +Prir 0, 9r)d, +P, +1(6,9,—5, 9,)4, (65) 
where g’, is the transform of g, under b, > b;. 
The scheme will be independent of sign if and only if 


(63) 


9,4, _ b,9,d, (66) 
and [e,{6,+ (1—2b,)(1—g;)}+-6,(9,—g;)]d, = 9. (67) 
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The possible solutions of (66) are 
(i) 9,=9 (ii) g,=1-—6b, (iii) g, = 1—d, (iv) g, =, 
g, = 0 g, = 140; g, = 1—d, (v) gi =b,. 
Substituting these solutions in turn in (67) it will be seen that (67) will 


only be satisfied by (i) and (iii) but not by (ii) or (iv). The former reduce 
(63) to the common form 


d.1 = (6,,,—6,)°(1 —Pp,4,)—P,.15,(1 —2b,.,)d,. (68) 
Since this is obtained by putting s,,,, = 0, (68) is the combined version 
of (37) and (41) and represents the method of short cut in the scale of four. 
Hence we have the result that this is the only scheme for use with 
multi-input type multipliers in a serial machine. 
In practice it may be more convenient to use the schemes (48) or (50) 
and precede the multiplication in cases of a negative multiplier with one 
subtraction of the multiplicand from the product register. 


11. The application of the mathematical analysis to division 
processes 

The usual restoring technique of division involves for n-digit numbers 
n trial subtractions and on average of }n corrections—a total of jn 
‘addition’ times. 

To reduce the division time von Neumann suggested that the non- 
restoring technique should be used when every division requires n addition 
times. This can be generalized to a process in which for each digital posi- 
tion the divisor is added to or subtracted from the remainder or not 
according to some rules. In this case the quotient is represented in the 
form ¥ (—)*d,2”" which has been studied in the earlier sections. 

Two questions then arise. Can the digits b be recovered from the form 
> (—)*d,2"? What is the rule deciding whether to change the remainder 
which will produce the minimum digit representation of the quotient? 

The first question is easily answered. Rearranging (14) we obtain 


d,,, = by, +61 —2b,,,)(1—d,) +8,d,(1—2b,,,) 
By ty sq(1—2by44) = b(1—d,) +8, 0, 
by .4(1—2d),4)-+dy,, = 6(1—d)) +8), 
bray = yxy +{b{1—d,) +-8,d,}(1 — 2d, ,,). (69)+ 


This defines the 6’s recursively from 6, = d,. In the special case of 
ordinary non-restoring division d, = 1 for all r, when (69) reduces to 


+ In the derivation of (69) free use is made of the identity (1—2z)* = 1. 
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This, with the forcing rule b, = 1 (for positive quotients) enables the 
digits of the quotient to be obtained in descending order as the division 
proceeds. This result was first obtained by von Neumann. 

For the genera! form, the recursion from most to least significant digits 
is not possible. Rearranging (69) we have 


b(1—d,) = d),,+b,.,(1—2d),,)—8)d). (71) 
If d, = 0 this gives by = dys +6,,,(1—2d,,,) 
= (141-41)? (72) 


but if d, = 1, (71) gives no information about 5,. In the division process, 
however, the relative signs of remainder and divisor determine the value 
of s for the next non-zero d and with this added information a reverse 
recursion is possible. 

Suppose, without loss of generality, positive dividend and divisor and 
record s = 0 for subtractions, s = 1 for additions. Let ¢,, denote the sign 
of the remainder at the end of a cycle to determine the digit d,, t = 0 
denoting positive remainder. Then 


ha tus (73) 
The fact that the sign of the remainder does not change unless an addition 
or subtraction is performed can be expressed by the equation 
(t,-1—t, )(1—d,,) = 0. (74) 
Consider the sequence 


d,, =1, ‘dy — dy», ooey dys =0, d,=1 


u > 
and let 6, = 1—zx. Then applying (49) for t = u—1 we obtain b,_, = 2. 
Similarly 6,_, = b,_3,..., 6,4, = x. Now apply (49) with t = wu; this gives 
X= By = bys = bysoreey ty = t, ie. b, = 1—t,. Similarly if 


d,, = dy4 =1, b,, - i~s,4 = 1-t 


u* 





Thus the recursive rule for obtaining the b’s is 
b, _ (d,.,;—6,.,)°(1—d,)+(1—t,)d,, (75) 
which reduces to (70) if d, = 1 for all r by using (73). 


We now turn to the second question. In the notation of equation (2) put 


u 


R, = | ¥ (—)»d, 2" |. 
' 0 


Then we have the result that the minimum digit representation is the only 
form which for all integers has the property: 


Rus = R, or Ry, >2R, for all u. (76) 
Clearly, if R.., * B,, Rua = Quel R.. 


u* 
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If the upper sign applies (76) clearly holds. If the lower sign applies, 
R,.,—2R, = 2%1—3R, 
= 2+1_ 3d, 2"43R, . (77) 
This must hold for all sequences b, and, in particular, for b, = 0 (r < u). 
For all s, this implies d, = 0 (r < uw) and hence R,_, = 0. In this case, 
d, = 0. Thus, d,d,,, = 0 is a necessary condition for (76). It is also 
sufficient, for in this case 


Ry, < x 4° = 4(4?+1—}), 


s=0 


R,, = }{2"+?—1—u(mod 2)}, 
R,.,—2R, > Ry.y—3R,-, > 14+(u—1)(mod 2) > 0. 
The condition d,,,d,, = © gives from (14) 
{(by1—by)*—(6y—8,)(1—20y )}dy = 0, 
{Ousa+8y(1—2by44)} = (1—2by 4s )(8y—busr)dy = 9, 
i.e. if é,, = 1, &, = 6... 
Hence (76) implies d’s given by (15) and (16), the minimal digital repre- 
sentation. 
If the dividend is s,, the divisor D, and the quotient is N given by (2) 


then &% = ND 


and the remainder at the u’(n—u)’ stage is 


n u 
8, = 8o—( | (—)*d,2")D = (> (—)*d, 2")D, 

s,| = R,|D}. (78) 

Hence for minimal digit representation 
Su4a| = |8y| OF [843] > 2\8,!. 
Thus the remainder changes if and only if 
is,,| > |D)2". (79) 

Thus in a practical divisor, after an addition or subtraction, the remainder 
should be shifted forward until it first exceeds 3 of the divisor. This is 
not a practical possibility for a parallel machine, but this rule is easily 
mechanized into a serial divisor which determines two or three digits of 
the divisor in each cycle. 

To obtain two digits in each cycle an adder is provided which produces 
3s,, in each cycle. This is fed to adder-subtractor units which produce the 
correctly signed versions of 3|s,,|—|D| and 3\s,,|—2\D)|. Two other adder- 
subtractors produce the correctly signed versions of |s,/—|D| and 
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\s,,|—4|D| which are fed into delay lines. s, is also fed into a delay line. 

At the conclusion of the cycle the sign digits of the outputs are inspected 
and the appropriate gate opened to feed the correct value of s,, into 
the trebling adder during the next cycle. These same sign digits determine 
the value of d,, and d,_, from which the values of b,, 6,_, can be found 
using (75). 

This mechanization depends on the fact that both possible subtractions 
will not take place in any cycle. If d,_, = 1, it is known that d,_, = 0 
and the delay lines can be adjusted to pass s,,_, to the trebling adder. 

To obtain three digits in each cycle additional adders must be incor- 
porated to form the correctly signed versions of 


3\s,|—24|D|, 3\s,|—14|D|, |s,|—}|D|, and |s,|—14|D| 
in addition to the others, so that the possibility of the digital pattern 
d, = 1, d,_, = 0, d,_, = 1 can be dealt with. 

For parallel machines, the most natural rule is to change the remainder 
if |s,,| > « independent of D. If « is chosen as 2“-! this corresponds to 
shifting the remainder to standard form after each change. 

An analysis of the resulting digital pattern is very complex and the 
author has been unable to solve the problem. It is conjectured that the 
average number of digits in the resulting representation of a uniformly 
distributed quotient is half the number of places, but this is based on 
a dubious hypothesis concerning the distribution of successive remainders. 


Conclusions 


The application of Boolean algebra to the problems of making fast 
multipliers and dividers in a binary computing machine results in several 
practical procedures which have hitherto not been explored. 
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